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NEXT-GENERATION EARTH RADIATION BUDGET INSTRUMENT 

CONCEPTS 


Katherine L. Coffey 
(ABSTRACT) 

The current effort addresses two issues important to the research conducted by the 
Thermal Radiation Group at Virginia Tech. The first research topic involves the 
development of a method which can properly model the diffraction of radiation as it 
enters an instrument aperture. The second topic involves the study of a potential next- 
generation space-borne radiometric instrument concept. 

Presented are multiple modeling efforts to describe the diffraction of monochromatic 
radiant energy passing through an aperture for use in the Monte-Carlo ray-trace 
environment. Described in detail is a deterministic model based upon Heisenberg’s 
uncertainty principle and the particle theory of light. This method is applicable to either 
Fraunhofer or Fresnel diffraction situations, but is incapable of predicting the secondary 
fringes in a diffraction pattern. Also presented is a second diffraction model, based on 
the Huygens-Fresnel principle with a correcting obliquity factor. This model is useful for 
predicting Fraunhofer diffraction, and can predict the secondary fringes because it keeps 
track of phase. 

NASA is planning for the next-generation of instruments to follow CERES (Clouds and 
the Earth’s Radiant Energy System), an instrument which measures components of the 
Earth’s radiant energy budget in three spectral bands. A potential next-generation 
concept involves modification of the current CERES instrument to measure in a larger 
number of wavelength bands. This increased spectral partitioning would be achieved by 
the addition of filters and detectors to the current CERES geometry. The capacity of the 
CERES telescope to serve for this purpose is addressed in this thesis. 
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1.0 INTRODUCTION 


The Earth/atmosphere system is evolving due to human activities and to natural events. 
The IPCC (Intergovernmental Panel on Climate Change) has concluded that “the balance 
of evidence suggests a discernible human influence on global climate through emissions 
of carbon dioxide and other greenhouse gases” [IPCC, 1997]. Activities such as 
deforestation, construction, biomass burning, and agricultural and industrial activities, as 
well as volcanic eruptions, all alter the composition of the Earth/atmosphere system by 
changing the planetary reflectance, adding aerosols to the atmosphere, and increasing and 
altering atmospheric gases. The Earth’s environment has been polluted as never before in 
the past century, and various independent measurements indicate that the Earth’s 
temperatures are changing, although how much is not exactly known. Burroughs [1997] 
reports that there has been a general warming trend of between 0.3 and 0.6 0 C in the 
Earth’s meteorological temperature throughout the twentieth century; most of this change 
concentrated in the period between 1920 to 1940 and since the mid- 1 970 ’s. 
Discrepancies exist between satellite and terrestrial measurements of Earth’s 
temperatures, and it is always possible that some of the detected changes may be 
attributed to the advances in measurement devices rather than truly being due to actual 
changes in the Earth’s temperatures. Finally, the possibility exists that observed changes 



Katherine L. Coffey 


Chapter 1.0 Introduction 


2 


are a natural part of the behavior of the dynamic Earth/atmosphere system. We approach 
the end of the twentieth century with more questions than answers with regard to the 
health of our planet. Scientists continue in the search for answers to the myriad questions 
about some of the most challenging issues that have ever faced our planet. 

There is no question that human activity has altered the composition of the Earth’s 
atmosphere. Human activity has sent huge quantities of pollutants such as carbon 
monoxide, sulfur dioxide, nitrogen oxides, hydrocarbons, and particulate matter into the 
atmosphere. Scientists have observed many notable changes in atmospheric gases in the 
second half of this century. Stratospheric ozone, important for its ability to absorb 
ultraviolet solar energy thus shielding the Earth, is being depleted due to human activity 
such as the emission of nitrous oxide and chlorofluororcarbons. One study conducted at 
Mauna Loa Observatory in Hawaii showed that the concentration of carbon dioxide in the 
atmosphere has increased by about 10 percent since 1958. Using the evidence of air 
bubbles trapped in polar ice, and recent observations such as those at Mauna Loa, some 
scientists estimate that the atmospheric concentration of carbon dioxide has increased by 
up to 25 percent since the early 1800’s. Although there are other sources of atmospheric 
carbon dioxide (i.e. the decay of vegetation and volcanic eruptions), it is believed that the 
major culprit is the burning of fossil fuels [Ahrens, 1992]. These are but a few of the 
various upsets to the balance of the Earth/atmosphere system brought on by human 
activity. 

Other evidence found in ice cores, as described by Burroughs, sways opinion in an 
opposite direction and adds weight to the intriguing question of how much detected 
changes are truly a result of human activity. These ice cores, taken from the Greenland 
ice sheet, reveal a history of huge fluctuations in the Earth’s climate that occurred 
independent of man’s influence. The ice cores show that the climate has been fairly 
stable over the past 1 0,000 years, but during previous years, as far back as 1 00,000 years 
ago, a picture of a highly erratic climate emerges. Assuming that the detected upward 
trend of less than one degree Celsius measured over this past century is an accurate 
figure, the earlier changes revealed by the ice cores were five to ten times greater, and 
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occurred over only a few years rather than over the span of a century. The evidence 
revealed by these Greenland ice cores has provided a new perspective in climatic 
thinking. Where it was previously believed that big changes in climate could only occur 
over a long period of time due to the thermal inertia of the oceans, it is now conceivable 
that the climate could undergo huge shifts in short periods of time. Because of such 
recent discoveries, some climatologists believe that these natural fluctuations must be 
strongly considered in the study of global warming. As emphasized by Burroughs, 
without a better understanding of the extent of past natural climatic change it is not 
realistic to plan on the basis that current changes are the consequence of human activities. 

The detected slight rise in the Earth’s temperature, the fact that the activities of mankind 
are altering the Earth’s environment as never before, and the knowledge that the Earth’s 
climate has undergone dramatic changes in the past independent of man’s influence, all 
combine to pose a unique challenge to the science community in determining whether the 
current warming is linked to human activities, or is the result of natural mechanisms in 
the climate. In the quest for answers to these questions, some of the most useful tools 
are GCMs, or Global Circulation Models, which can be used to study cause-and-effect 
scenarios in the Earth/atmosphere system [Wielicki, et al., 1995]. GCMs are computer 
models that mathematically model the extremely complex physics of the 
Earth/atmosphere system. These models incorporate many approximations and 
simplifications, and there is much room for improvement as many unknowns still exist. 
The importance of GCMs which can yield trustworthy conclusions as to the cause/effect 
of global warming is paramount. Information revealed by these GCMs is used by 
organizations such as the IPCC (Intergovernmental Panel on Climate Change), which 
serves to assess scientific information about climate change relevant for the formulation 
of international and national policy. Because of their importance, there is much interest 
in the refinement of GCMs. Refinement can be accomplished by utilizing findings from 
continued studies of the Earth/atmosphere system. 

National and international efforts are being made to address these uncertainties. NASA’s 
most recent contribution to these efforts is the Mission to Planet Earth, which is part of 
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the United States Global Change Research Program (USGCRP). This research effort 
consists of a series of space-based remote sensing platforms, the largest of which is the 
Earth Observing System (EOS) [Anon, 1993]. EOS involves a series of Earth-orbiting 
satellites containing a variety of instruments designed to provide critical global 
observations of the Earth/atmosphere system, computing facilities, and science 
researchers, and is committed to data collection for at least a 15-year period. The EOS 
system is particularly effective, because a single satellite contains multiple instruments 
measuring many independent physical processes that can be observed simultaneously for 
a given scene type. Chapter 2.0 describes several NASA instruments, the measurements 
that they are used to obtain, and the data products they provide. This will serve as 
background to describe the current research effort. Figure 1.1 summarizes the many 
interrelated factors involved in the study of global wanning and climate change. 
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2.0 EFFORTS TO MONITOR GLOBAL CLIMATE CHANGE 

2.1 The Earth radiation budget 

Averaged over the entire Earth and over the span of a year, the sun delivers 
approximately 340 W/m 2 of shortwave radiation. Approximately 99 percent of this 
energy is below 3.8 pm, with one percent below 0.3 pm. Of the arriving solar energy, 
about 30 percent, or 100 W/m 2 , is reflected back into space. This reflected energy is 
called the Earth’s albedo. The solar energy which is not reflected back into space (240 
W/m 2 ) is absorbed by the Earth/atmosphere, which heats up and emits longwave 
radiation. The Earth’s temperatures currently range between approximately 240-300 K 
with 99 percent of its emitted energy at wavelengths longer than 4.6 pm and only two 
percent above 60-70 pm [Lenoble, 1993]. Ahrens [1992] explains that “Earth’s 
temperature” can mean a variety of things. For instance, the Earth’s observed average 
surface temperature is about 288 K. On the other hand, if Earth is viewed as behaving as 
a blackbody, where it absorbs solar radiation and emits infrared radiation at equal rates, 
its radiative equilibrium temperature is approximately 255 K.. Thus care should be taken 
when reporting and interpreting changes in the “Earth’s temperature”. 


6 
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Since the solar and terrestrial radiation are present in wavelength bands which barely 
overlap, the two can be measured separately. By subtracting the total emitted and total 
reflected energy from the solar contribution, the net radiation budget, or the net amount 
absorbed by the Earth at a given time and location can be computed. If the net amount at 
a given time and location is positive that particular area of the Earth tends to warm, 
whereas if the net is negative that area cools. This spatial variation in temperature is 
what drives the weather of our planet. If this difference is averaged temporally and 
spatially, the net radiation budget is approximately zero. However, if some disturbance 
(due to human activity or to natural events) were to occur which disturbed this 
equilibrium, a time-dependent climate change would occur, and the temperature of the 
Earth would change (increase or decrease) until radiative balance was once again 
established. One of the goals of Earth-observing instruments is to determine whether this 
sort of imbalance is occurring, and if so, to identify its causes. One of the steps in this 
process is to measure the parameters necessary to derive the global net radiation energy 
budget. 


2.2 Earth-observing instruments 

The use of satellites is an especially effective measurement method in the quest for the 
answers to the fundamental question of whether the temperatures of the Earth are 
changing, as they can provide global coverage and thus make measurements used to 
derive the global TOA (Top-of-the-Atmosphere) radiant energy budget. Satellite-based 
observation began in the early 1960’s with the onset of the Space Age, and has continued 
ever since. 

As the generations of Earth-observing instruments have developed, they have served in 
the continued gathering of data initiated by previous instruments, in making increasingly 
advanced measurements and in using increasingly sophisticated algorithms in data 
reduction. As with most improvements there are trade-offs, as these upgrades do not 
come without a price. In the process of replacing one instrument with another having 
greater capability, the organization leading the investigation risks the disturbance of the 
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collection of continuous science data. Detected changes in measurements which may be 
due to the transition in instruments could be misinterpreted as changes in the Earth’s 
energy budget. This risk is minimized as much as possible, and instruments are replaced 
as the motivations for improvement outweigh the risks. To help minimize risk, new- 
generation and past-generation instruments’ missions are made to overlap so that 
calibrations can be transferred, if at all possible. Instrument improvement is achieved by 
using the lessons learned from previous missions, and the progressive state of knowledge 
about the Earth/atmosphere system that reveals data products which may be more useful. 
Change is also motivated by the organization’s desire to stay atop the latest technology. 
In the case of the National Aeronautics and Space Administration, change is often driven 
by the fact that “NASA Headquarters and Administration appear to want to view NASA 
as technology developers rather than scientific leaders” [Barkstrom, 1998]. So for 
various reasons, the push is on to conduct successful experiments using state-of-the-art 
technology and to gather useful information that can be used to advance our knowledge 
about the Earth/atmosphere system while planning constantly for the future generations 
of instruments to be used to monitor our planet. At the same time, Earth-observing 
experiments must continue to gather a complementary set of data throughout the 
generations of instruments so that the studies are sufficiently long-term to monitor change 
over the decades. 


2.3 Earth-observing instruments leading to the current research effort 
2.3.1 ERBE 

In 1979, the National Aeronautics and Space Administration began ERBE (the Earth 
Radiation Budget Experiment), a mission that was to measure the most basic parameters 
in monitoring the Earth’s climate. ERBE was launched aboard the NASA ERBS (Earth 
Radiation Budget Satellite) in 1984, and aboard NOAA-9 and NOAA-IO in 1984 and 
1986. The ERBE instrument is a scanning thermistor bolometer radiometer consisting of 
three channels sensitive in the short (0.2-5. 0 pm), the long (5.0-50 pm) and the total (0.2- 
100 pm) wavelength bands. The principal goals of this mission were to measure 
broadband radiances at the Top-of-the-Atmosphere (TOA), to convert these anisotropic 
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radiances to TOA fluxes using Angular Distribution Models (ADMs) and the Maximum 
Likelihood Estimation (MLE) technique [Wielicki and Green, 1989], and to derive the 
global TOA radiation energy budget. ERBE collected thirteen years of data, and 
provided what NASA calls the most accurate data of the Earth’s outgoing longwave and 
solar reflected shortwave radiation ever obtained, as well as answers to some long- 
standing questions about climate forcing and feedback mechanisms in the 
Earth/atmosphere system. Among these important findings, ERBE data showed that the 
annual average effect of clouds is to cool the current climate system. This did not, 
however, end the debate as to whether clouds act to decrease global warming. “A 
common misconception is that because clouds cool the present climate, they will likewise 
act to moderate global warming. What is actually important is the change in the net 
cloud radiative forcing, associated with a change in climate, that governs cloud feedback” 
[Weilicki et ai, 1995]. In the process of answering questions, these answers prompted 
further questions and the next generation of instruments evolved. The insight into the 
influence of clouds gained by ERBE data served as a basis for naming the monitoring of 
clouds and their influence on the Earth’s radiant energy system as a top priority in the 
next generation of instruments, CERES. 

2.3.2 CERES 

The Clouds and the Earth’s Radiant Energy System (CERES) is a suite of broadband 
scanning radiometers based on the ERBE instrument, but featuring many improvements 
[Wielicki, 1996; Bongiovi, 1993; Haeffelin, 1996; Priestley, 1997; Smith, 1998], The 
CERES instrument incorporates full two-dimensional directional sampling by scanning in 
both elevation and in azimuth angle. It includes the same short and total wavelength 
channels of ERBE, but the longwave channel was replaced by a “window” channel, 
sensitive in the region of 8.0-12.0 pm. The first CERES instrument, the PFM (Proto- 
Flight Model), was launched in late 1997 aboard TRMM (Tropical Rainfall Measuring 
Mission), as part of NASA’s Earth Observing System (EOS). Two more CERES 
instruments are scheduled to be launched aboard the spacecraft EOS-AM and EOS-PM in 
the late 1990’s and the beginning of the 21 st century, and these and follow-on instruments 
will extend measurements for a total of fifteen years. CERES data will serve to extend 
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the thirteen years of ERBE data by measuring broadband radiances at the top of the 
atmosphere which will be converted to TOA fluxes. The CERES Pathfinder Project was 
organized to analyze CERES data products and to “bridge the gap” between ERBE and 
CERES missions. In addition to TOA fluxes, more sophisticated parameters are being 
determined. The CERES instruments fly on satellites containing cloud imagers that make 
simultaneous measurements of the same scene being viewed by CERES. CERES data 
are being used with this cloud imager data to infer surface fluxes and to provide the 
vertical profile of radiative divergence. Data collected during CERES azimuthal 
rotations is being used to build angular distribution models (ADMs), striving to meet the 
mission goals to reduce the ADM errors by a factor of four over ERBE. CERES is 
expected to provide TOA fluxes that are two to three times more accurate than those of 
ERBE data [Wielicki et al., 1996]. Over the course of the development and launch of 
CERES, and with the use of CERES Science data, NASA has identified the desirable 
features of the next generation of instruments. In summer, 1998, NASA outlined these 
features and named this potential future instrument, PERSEPHONE. 

2.3.3 PERSEPHONE 

NASA developed the preliminary conceptual design of PERSEPHONE under the 
criterion that the next-generation instrument be smaller, less resource intensive, less 
costly, and requires less build time. In the spirit of these requirements, and in the 
continued quest for the understanding of climate forcing and feedback mechanisms, 
NASA has identified increased spectral partitioning of radiance as the principal feature of 
this potential next-generation instrument. This partitioning will involve dividing the 
measured energy into a larger number of spectral bands rather than simply measuring in 
the shortwave, total, and window channels as in ERBE and CERES. Since radiant energy 
is spectrally and spatially altered as it passes through the Earth/atmosphere system before 
arriving at the detector, and since most atmospheric constituents and surface properties 
are selective absorbers (i.e. they strongly affect radiation in limited spectral regions), 
measurement of energy in these narrower spectral bands will yield a deeper 
understanding of our Earth’s climate system [Barkstrom, et al., 1998]. 
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In order to achieve the goal of refined spectral partitioning while keeping cost at a 
minimum, NASA has proposed that the current CERES telescope be modified to serve 
several detectors. Each detector would measure energy in a different spectral band by 
placing the desired filters in its optical path. While expanding the capability of the 
current CERES telescope, NASA wishes to improve the calibration accuracy and spatial 
and angular sampling capability while maintaining the current envelope of size, mass, 
and electrical power. In addition, NASA wishes to maintain the same (or better) quality 
of the Optical Point Spread Function (OPSF) on all detectors as that of the current 
CERES telescope. The OPSF describes the radiation throughput to the detector as a 
function of the angle at which collimated radiation enters into the instrument. These 
concepts will be described more thoroughly in Chapter 5.0. 

Modification of the CERES telescope to meet the criterion outlined above creates new 
challenges that will require the state-of-the-art in technology to overcome. Barkstrom, et 
al. [1998] identify the three most serious challenges that must be overcome in order for 
this redesign to succeed. The first issue involves the need for an increase of detector 
sensitivity by more than an order of magnitude over the current CERES detector since the 
redesign calls for the division of energy into smaller spectral intervals. The other two 
issues involve the need to isolate the different spectral bands without loss of calibration 
accuracy, and the need for calibration sources for the narrower spectral intervals. 

In addition to these issues, a redesign of the optics may be required in order to achieve an 
approximately uniform radiative flux over all detectors. One of the topics of research for 
the current master’s thesis involves the determination of whether a redesign is required, 
and if so, further study of a potential optical prescription that will yield acceptable 
performance. 


2.4 The Thermal Radiation Group 

The Thermal Radiation Group of the Mechanical Engineering Department at Virginia 
Tech has been involved in the numerical modeling of Earth-observing instruments since 
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the early 1970’s, in addition to other research projects. Much of this work was funded 
under NASA grants, in support of the work done by the Radiation Sciences Branch, 
Atmospheric Sciences Division of NASA Langley Research Center. Under the direction 
of Dr. J. R. Mahan, master’s and Ph.D. students of the Thermal Radiation Group have 
worked to develop high-level, dynamic, electrothermal, end-to-end numerical models of 
both ERBE and CERES instruments. These models are capable of simulating the 
response of these instruments to simulated Earth scenes [Haeffelin, et ai, 1997], The 
end-to-end models include an optical/thermal-radiative module and a thermistor 
bolometer dynamic electrothermal module. The optical/thermal-radiative module is used 
to model the optics of the instrument using the Monte-Carlo ray-trace (MCRT) method as 
it applies to radiation heat transfer. The output from this module is then used as input to 
the second electrothermal module, a finite-difference or finite-element code 
characterizing the electrothermal behavior of the detector to various inputs. 

In addition to this work, over the course of the past few years a Ph.D. student in the 
Thermal Radiation Group, Felix Nevarez, has been developing a flexible, ray-trace-based 
numerical tool that can be used to build the radiative model of any instrument. In the 
past, students have built a highly specific code for the instrument at hand, such as those 
described for the CERES and ERBE instruments. This new tool brings with it the 
capability to build in a matter of weeks what could have taken up to several years using 
previous capabilities. The optical prescription for the PERSEPHONE instrument has 
been studied using this tool. 

In addition to this topic of investigation, the author has investigated an entirely separate 
issue, the proper modeling of diffraction of radiant energy as it enters an instrument as 
treated in the Monte-Carlo ray-trace environment. This issue is important in many of the 
computer modeling efforts conducted by the Thermal Radiation Group, but has been 
neglected in the past. It is believed that a higher level of sophistication in the modeling 
of diffraction will benefit modeling efforts by future generations of students of the 
Thermal Radiation Group. 
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2.5 Goals of the current research 

The goals of the current research effort are to: 

1. Develop a radiative model of the CERES instrument using the new ray-trace tool 
being developed in the Thermal Radiation Group. This model will be used to: 

■ serve as validation of the new ray-trace tool by comparing results to results from a 
previous radiative model and experimental results. 

■ study the current CERES telescope to determine whether changes will be required 
for the successful implementation of the multiple detector concept for potential 
next-generation instruments. 

■ investigate a new optical prescription which may better serve the needs of the 
proposed next-generation instrument concepts. 

2. Investigate the importance of the modeling of diffraction of radiant energy as it enters 
instruments of the type modeled by the Thermal Radiation Group. 

3. Investigate, develop, and test means of numerically modeling diffraction in the 
Monte-Carlo ray-trace environment. 



3.0 PROBLEMS IN OPTICS 

Much of the work done by the Thermal Radiation Group involves the development of 
radiative models of radiometric instruments. The Monte-Carlo ray-trace (MCRT) 
method is generally used, but the design objectives vary with the measurements being 
sought. For instance, the design of a solar aureolemeter used to measure the 
concentration and size distribution of atmospheric aerosols is a continuing topic of 
investigation in the Thermal Radiation Group. This is an imaging instrument that 
measures the circumsolar radiation pattern produced by atmospheric aerosols. These 
measurements are important because the scattering produced by aerosols in the region 
around the sun is much greater than that due to other atmospheric constituents. The 
image quality is of principle concern in the design of this instrument [Deepak and Adams 
1983; Nakajima, et al., 1983], On the other hand, the CERES instrument is not an 
imaging instrument, but a radiometer. Energy throughput, rather than the image quality, 
is the principle concern in the design of this instrument. As described by Walkup [1996], 
the area in which the two disciplines of imaging and radiometry merge is imaging 
radiometry, illustrated in Figure 3.1. This best describes the type of investigations 
conducted within the Thermal Radiation Group. 


14 
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Figure 3.1 Combination of imaging and radiometry disciplines, (borrowed from Walkup 
[1996]). 


3.1 Ideal versus real behavior in optical systems 

Much of the work presented in this thesis in some way pertains to the deviation of the 
true performance of instruments from simplified models of their ideal behavior. The 
ideal behavior of optical systems is often called Gaussian optics. This first-order theory 
is used for preliminary design specifications, and can be used for mirror placement and 
for locating the ideal image point in an optical system. Gaussian optics describes a point- 
to-point, object-to-image relationship. That is, when rays leave a point within the object 
plane, they will arrive at a single point in the image plane. True optical performance of 
instruments deviates from this idealized behavior. The degree of this deviation can be 
described by five different aberrations: spherical aberration, coma, astigmatism, field 
curvature, and distortion. For the purposes of the current effort, it is important to 
describe spherical aberration in an optical system. For a detailed description of the other 
aberrations, refer to Walkup [1996]. 

3.1.1 Spherical aberration 

Spherical aberration is the only aberration that occurs on-axis. Because of this 
aberration, rays reflected from conic reflectors do not converge to the ideal infinitely 



Katherine L. Coffey 


Chapter 3.0 Problems in optics 


16 


small spot size as they cross the optical axis at the Gaussian image plane. Instead, a blur 
circle of finite size is formed at the focal plane. Rays that traverse the telescope near the 
optical axis behave ideally, passing through the Gaussian image point, while those 
reflecting far from the optical axis do not. The point at which this blur circle is a 
minimum is called the “circle of least confusion,” as illustrated in Figure 3.2. It is 
important in the design and assembly of an imaging telescope that this circle of least 
confusion fall in the detector plane. For a non-imaging telescope such as CERES, the 
telescope is assembled such that the circle of least confusion is located at the plane 
containing the precision aperture that serves to shape the image cast on the detector 
beneath. 



Figure 3.2 Illustration of spherical aberration (borrowed from Walkup [1996]). 
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3.2 Ideal versus real propagation of electromagnetic waves 

The idealized behavior of the propagation of electromagnetic energy is sometimes 
referred to as the ideal ray approximation. This model ignores diffraction, and rays are 
modeled as straight lines that are perpendicular to the propagating wavefronts. As 
electromagnetic waves enter apertures or pass by obstacles, the ideal ray approximation 
does not account for diffraction effects that are known to occur. The correction of this 
idealization in instrument modeling is a major topic of this thesis. The following 
discussion of diffraction serves to provide the necessary background for understanding 
the work presented in Chapter 4.0. 

3.2.1 Diffraction 

3.2.1. 1 Basic ideas 

Diffraction and interference are actually one and the same: the constructive and 
destructive interference of in-phase or out-of-phase waves arriving at the same point. 
The convention has developed to call this additive wave property interference when it 
deals with only a few waves, and to call it diffraction when it involves many waves 
[Hecht and Zajac, 1974], Diffraction occurs when electromagnetic waves pass through 
apertures or past obstacles, causing the waves to deviate from their original lines of 
travel. The degree to which this divergence occurs depends upon the wavelength of the 
entering radiation relative to the aperture dimensions. For the case where the wavelength 
of the approaching wave is much less than the aperture dimensions, diffraction effects are 
slight. In this case, the ideal ray approximation is valid. The resulting intensity on an 
observation screen placed behind the aperture is in the shape of a top hat, with no light 
arriving outside the area which is formed by the projection of the aperture area on the 
screen. This case is illustrated in Figure 3.3 (a). When the wavelength of the light is on 
the order of the aperture size, diffraction effects are more pronounced as shown in Figure 

3.3 (b). If the entering energy is monochromatic, a pattern of alternating minima and 
maxima will form on an observation screen placed behind the aperture. Energy will 
spread out to areas that would be predicted to be in shadow by the ideal ray 
approximation. The shape and span of the resulting diffraction pattern depends on the 
aperture size, the wavelength of the entering energy, and the distance between the 
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observation screen and the aperture. Finally, when the wavelength of the approaching 
energy is much greater than the aperture dimensions, the aperture behaves as a point 
source that emits spherical waves, (i.e. a slit behaves as a line of point sources) as 
illustrated in Figure 3.3 (c). In this case, the energy will spread out in the region where a 
shadow would be predicted, but no pattern of maxima and minima will be observed. 



Figure 3.3 The behavior of radiation as it passes through an infinite slit: (a) the 
wavelength of entering radiation is much less than the slit width, (b) the 
wavelength is approximately equal to the slit width, and (c) the 
wavelength is much greater than the slit width (borrowed from Serway 
[1994]). 


The diffraction models that are presented in Chapter 4.0 are used to simulate the behavior 
of monochromatic, coherent radiant energy as it passes through an aperture, and the 
results are compared to closed-form analytical solutions of the expected interference 
pattern for validation. Upon validation, these models will then be used to simulate 
diffraction for cases in which the approaching energy is neither coherent nor 
monochromatic. Such is the case for applications such as instruments used to measure 
the Earth’s radiation energy budget. 


3.2.1.2 The two diffraction regimes 

Diffraction can be classified into two regimes: Fresnel, or near-field, and Fraunhofer, or 
far-field, diffraction. Fraunhofer diffraction occurs when the light source and the 
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observation screen upon which the diffraction pattern is observed are both effectively at 
infinite distances from the aperture causing the diffraction. Fresnel diffraction occurs 
when either (or both) the source or the screen are at finite distances from the aperture. In 
the development of a radiative model of an instrument, one must determine how crucial 
the modeling of diffraction is. The first step in this assessment is to determine which of 
the two diffraction regimes governs the diffraction. For the case in which energy is 
entering through a rectangular slit, and the source is at infinity, the parameter 

i 

{ 7 \ 2 

A^ = af 2 , (3-1) 

\KZJ 

where a is the slit width, z is the distance between the aperture and the observation 
screen, and k is the wavelength of the entering energy, is used for the purpose of 
categorizing diffraction as Fraunhofer or Fresnel [Wyatt, 1987; Haskell, 1970]. At, =1 is 
taken to be the transitional configuration between the two regimes. For cases in which 
A^>1, Fresnel diffraction prevails, while for cases in which A£<1, Fraunhofer diffraction 
prevails. Figure 3.4 illustrates this transition, brought on by moving the observation 
screen away from the aperture. In this case, the point source is at infinity with respect to 
the slit. As the observation screen is moved away from the aperture, the rays become 
approximately parallel upon their arrival at the observation screen, meeting the criterion 
for Fraunhofer diffraction. Note that a similar transition would occur if the wavelength 
of the entering radiation were increased while holding the parameters z and a constant. 
Such is the case in the practical application presented in Chapter 4.0. 


The diffraction patterns shown in Figure 3.4 are based on the closed- form analytical 
solution reported by Haskell [1970]. The intensity along the screen is given in terms of 
the Fresnel integrals 


and 



(3.2 a) 



(3.2 b) 
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and the parameters 



(3.3) 

5> = ia4[(5./A«)+*]. 

(3.4) 

where is as defined in equation 3.1 and 


^ = - 2x lrrl 2 - 

(3.5) 


Thus the intensity along an observation screen, as illustrated in Figure 3.4 is given by 


I(x) = i [C(l;0 - C«of + i[S(50 - SffiOp . 


(3.6) 


Equation 3.6 can be used to determine diffraction patterns in either the Fresnel or 
Fraunhofer regime. However, Fraunhofer diffraction can be described by simplified 
equations found in Laufer [1996]. 
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Monochromatic, coherent radiation of 
wavelength, X = 0.001 pm approaching 
infinite slit of width, a = 0. 1 pm from a 
source at infinity. 

X 





A$=1.0 

x, distance from central axis, pm 





Figure 3.4 Illustration of the transition from Fresnel to Fraunhofer diffraction and the 
resulting change in the intensity distribution, I(x) where y = I(x) 
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Table 3.1 summarizes the characteristics of the two diffraction regimes. The modeling 
capabilities mentioned refer to the models that are presented in Chapter 4.0. 


Table 3.1 Comparing characteristics of Fraunhofer and Fresnel diffraction. 


Fraunhofer 

Fresnel 

♦ Far-field 

♦ Near-field 

♦ Simpler analytical solutions. 

♦ Complex analytical solutions. 

♦ Instrument throughput not good. 

♦ Instrument throughput good. 

♦ Possibility of modeling phase using new 
ray-trace approach (Model 2) and MCRT 
method. 

♦ Cannot model phase using MCRT 
method. 

♦ Good modeling is more important, ideal 
ray trace deviates greatly from , actual 
behavior (Figure 3.5 (a)). 

♦ Good modeling is not always as 
important, ideal ray trace approximates 
the true behavior well in very near field 
(Figure 3.5 (b)). 



(a) (b) 

Figure 3.5 Illustration of the deviation of the ideal ray approximation from the true 
intensity pattern resulting from the passage of energy through a slit 0. 1 pm 
wide (a) for a Fraunhofer diffraction situation, and (b) for a Fresnel 
diffraction situation. 

Figures 3.5 (a) and (b) illustrate the importance of modeling diffraction well in the 
Fraunhofer regime, and the reduced importance of modeling diffraction well in the 
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Fresnel regime. Because energy is distributed far into the geometrical shadow in 
Fraunhofer cases, much of the on-axis energy is dispersed out of the optical path. For 
this reason, Fraunhofer diffraction is generally avoided in imaging systems. Ideally, an 
effort should be made to ensure that diffraction occurs in the Fresnel region in the early 
design stages of an instrument [Wyatt, 1987]. 

3.2.1.3 A brief history of diffraction 

The diffraction of electromagnetic waves has puzzled and challenged scientists, engineers 
and physicists for more than five centuries, and some of these challenges remain even 
today. The authors of most optics and physics books state that diffraction is one of the 
most complicated problems to be dealt with in optics. Diffraction is of academic 
importance, as it serves as evidence of the wave nature of light. It is of practical 
significance to engineers and scientists as it limits the resolving power of optical devices, 
and often results in less than ideal instrument throughput. 

Challenges posed by the phenomenon of diffraction are interwoven throughout the 
history of optical science. Every generation of scientists since the 1400s has made a 
contribution to further our understanding of diffraction. The following briefly 
summarizes the stages of diffraction knowledge, borrowing from the historical summary 
by Hecht and Zajac [1974], and other authors. Throughout this historical summary, many 
important principles are described which serve as background for the diffraction models 
presented in Chapter 4.0. 

The first reference to the phenomena of diffraction was made in the work of Leonardo Da 
Vinci (1452-1519). It was not until 1665 that diffraction was accurately described by 
Professor Francesco Maria Grimaldi in his observations of bands of light within the 
shadow of a rod illuminated by a small source [Bom, 1970]. Soon after, Robert Hooke 
(1635-1703) performed experiments in which he also observed diffraction patterns. 
These observations of diffraction prompted the beginning of the wave theory of light, and 
the great debate as to whether light is corpuscular (stream of particles) or wave-like. 
Christian Huygens (1629-1695), a proponent of the wave theory, continued to extend the 
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wave theory of light. One of his significant contributions, Huygens’ principle which 
states that “every point on a primary wavefront serves as the source of spherical 
secondary wavelets such that the primary wavefront at some later time is the envelope of 
these wavelets” is illustrated in Figure 3.6 (a). This model, which treats the secondary 
wavelets as radiating equally in all directions, and the envelope of the wavelets as the 
intensity at an observation screen behind the aperture, is unable to predict the expected 
interference pattern due to the diffraction of energy. It also predicts a backward 
propagating wave, which is not observed in reality. Thomas Young (1773-1829) was the 
next important contributor to the wave theory of light despite ridicule from many of his 
colleagues who supported the corpuscular theory. Young provided the first clear 
demonstration of the wave nature of light in 1801. His double slit experiment showed 
that light coming from a single source, arriving at a point by two separate paths, 
combines both constructively and destructively, resulting in an interference pattern on an 
observation screen placed behind the two slits, as illustrated in Figure 3.6 (b). 



Figure 3.6 (a) Illustration of Huygens’ principle for plane waves propagating to the 

right (excluding the backwards-propagating wave), and (b) illustration 
of Young’s double-slit experiment. 


Augustin Jean Fresnel (1788-1827) added to Huygens’ ideas in order to provide a model 
which could predict the intensity variations expected by diffracted energy. He added an 
undefined angular dependence factor, K, called the obliquity factor, which serves to 
model the directionality of the emission of secondary wavelets coming from a point 
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source. The Huygens-Fresnel principle states that “every unobstructed point of a 
wavefront, at a given instant in time, serves as a source of secondary wavelets of the 
same frequency as the primary wave. The amplitude of the optical field at any point 
beyond is the superposition of all of these wavelets (considering their amplitudes and 
relative phases)”. Fresnel began to describe diffraction on a mathematical basis and was 
able to predict diffraction patterns of light when it passed by obstacles. Kirchoff was the 
first to put the conceptual Huygens-Fresnel principle on a more mathematical basis, 
proposing a specific form of the obliquity factor, K which is used in Chapter 4.0. Instead 
of conceptual arguments to describe diffraction phenomena, mathematical descriptions 
began to be sought and used more extensively. The determination of an exact solution 
for a particular diffracting configuration is amongst the most difficult to be dealt with in 
optics [Hecht and Zajac, 1974]. The first such solution was published in 1896 by Arnold 
Sommerfeld (1868-1951). This solution utilized the electromagnetic theory of light, and 
involved an infinitely thin, opaque, perfectly conducting plane screen. Note that rigorous 
solutions such as the one published by Sommerfeld exist for very few configurations of 
practical interest today. Because of this, numerical methods based on conceptual 
diffraction models and other approximations are often used to determine the effect of 
diffraction on configurations of interest. 

3.2. 1.4 Modern methods of dealing with diffraction 

A number of methods for modeling and/or characterizing diffraction in an optical 
instrument have been developed in this century, taking advantage of advancements in 
computing power. Some of these methods are described briefly in the following 
paragraphs, and references are cited which provide more extensive details. 

Airy disk construction 

Circular baffles, apertures, or obscurations are quite common in optical systems. For 
such configurations, the importance of diffraction in their optical performance can be 
approximated by a simple calculation. The size of the diffraction blur or Airy disk, 
named in honor of British mathematician, Lord George Biddell Airy (1801-1892), is 
given by 
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Di = 2.44 X f , (3.7) 

where Di is the diameter of the diffraction blur, X is the wavelength of the entering 
energy, and f is the focal length of the optical system divided by the aperture diameter 
[Riedl, 1997], For the case of energy entering through a circular aperture, Di defines the 
area on the observation screen containing 84 percent of the incident energy, assuming no 
other aberrations are present. The diameter, D 2 , which defines the area that will contain 
91 percent of the incident energy is given by 

D 2 = 4.48 X f . (3.8) 

It is important to note in instrument design that there is no point in decreasing the 
diffraction blur size below the size of the blur circle due to the geometrical aberrations 
present in the optical system. Likewise, an optical system can never perform better than 
the limits that are predicted by diffraction, thus optical designers may reach a point where 
further minimization of geometrical aberrations would be in vain. Another important 
situation, a circular aperture with a central obscuration, is commonly found in telescopes 
such as CERES. Details of the treatment of this situation are provided by Riedl [1997] 
and by Walker [1998]. Note that the addition of a circular obscuration in the center of a 
circular aperture causes the incident energy to become more dispersed, thus the Airy disk 
diameters are larger. The degree of this dispersion of energy out of the desired optical 
path depends on the diameter of the obscuration relative to the diameter of the aperture, 
characterized by the percent obscuration. The larger the percent obscuration, the more 
the incident energy is dispersed. For this reason, systems with obscurations greater than 
50 percent are rarely used [Walker, 1998]. For further discussions regarding the 
influence of diffraction on instrument imaging performance, the interested reader is 
referred to Braun [1970], Tschunko and Sheehan [1971], and Boivin [1977], 

Method of moments 

The method of moments can be used to approximate the diffraction pattern of energy 
entering an aperture. This method, based on the method of weighted residuals, is a 
common tool within the electrical engineering culture, and can be applied to treat the 
propagation of electromagnetic waves of any wavelength. For details of this method see 
Hubing [1991], Graves et al. [1976] address the modeling of the electromagnetic field 
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penetration through a circular aperture in both the near and far field using this method. 
Results are shown to compare well with theory. 

Geometrical theory of diffraction 

The geometrical theory of diffraction provides a means of supplementing standard 
geometric optics with the ability to model diffraction. Keller [1957, 1962] and Balanis 
and Peters [1968] published descriptions and sample applications of this method, which 
treats diffraction as an edge effect, as originally proposed by Thomas Young [Keller, 
1962]. In addition to the usual rays, the geometric theory of diffraction introduces new 
diffracted rays which originate as a ray strikes an edge or vertex. As an example of the 
application of this method, the treatment of radiation as it enters an aperture is described 
by Keller [1957, 1962]. Rays which strike an edge become diffracted rays. When the 
incident ray is normal to the edge, the* diffracted rays produced are also normal to the 
edge (i.e. all lie in the same plane), and leave in all directions within that plane, as shown 
in Figure 3.7 (a). When waves approach the edge obliquely, the diffracted wave fronts 
are defined by cones emanating from the edge, as illustrated in Figure 3.7 (b). These 
obliquely diffracted rays follow the law of edge diffraction , which states that a diffracted 
ray and the corresponding incident ray make equal angles with the edge at the point of 
diffraction. The incident and diffracted rays lie on opposite sides of the plane normal to 
the edge at the point of diffraction [Keller, 1962]. All rays are assigned an initial phase 
upon their arrival to the aperture opening, and their final phase upon arrival to an 
observation screen is proportional to the optical length of a ray from this point. Keller 
[1957, 1962] and Balanis and Peters [1968] demonstrate good agreement when 
comparing results from the application of the geometrical theory of diffraction with the 
closed-form analytical solution of the expected diffraction pattern for several cases. 
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Knife edge 


Incident ray 



Figure 3.7 Illustration of the formation of diffracted rays, (a) when an incident ray 
strikes perpendicular to the edge, and (b) when incident rays strike the 
edge obliquely (borrowed from Keller [1962]). 


The geometrical theory of diffraction was implemented in the GUERAP II code, a 
computer program for the analysis of the stray radiation rejection capabilities of optical 
systems, as briefly described by Boyce [1977]. A subsequent version of this code, 
GUERAP III, replaced the use of this methodology with a newer technique, referred to in 
this thesis as the statistical method for dealing with diffraction [Morbey and Hutchings, 
1993]. 


Statistical method for dealing with diffraction 

This method is based upon Heisenberg’s uncertainty principle and the particle theory of 
light, and is convenient for use in the Monte-Carlo ray-trace environment [Chou, 1974], 
This approach does not keep track of phase, so the resulting pattern approximates the 
minima and maxima which occur alongside the central maxima. Heinsch and Chou 
[1971] and Likeness [1977] briefly describe the concepts of this approach, omitting most 
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of the details necessary for its implementation. The details needed for the application 
have been rediscovered by the author of this thesis and are discussed at length in Chapter 
4.0. 

Use of Huygens’ principle to model diffraction 

The use of Huygens’ principle to model an aperture as many point sources emitting 
isotropically in all directions has been implied by a technique described by Sinnott 
[1987], The diffraction pattern resulting from the passage of starlight through a telescope 
is desired. The Monte-Carlo technique is applied, sending random parallel rays through 
the telescope. Ray path lengths are determined upon their arrival at an observation 
screen, and the net phase at a given point is given by the sum of the phase of all rays 
arriving at that point. Results are given as a two-dimensional scattergram, and are shown 
to compare qualitatively with analytical results (i.e. the spatial distribution of energy 
compares with the expected diffraction distribution). Results are not, however, 
quantitatively compared with theory (i.e. comparison of the relative magnitude of 
secondary maxima to central maxima is not made). Further discussion and modification 
of this technique are proposed in Chapter 4.0. 



4.0 MODELING DIFFRACTION IN THE MONTE-CARLO RAY- 
TRACE ENVIRONMENT - 

In the Theory of Heat Radiation. Max Plank [1959] states that “so far as their physical 
properties are concerned, heat rays are identical with light rays of the same wavelength”. 
His book, which contains much of the basis upon which modem radiation heat transfer is 
based, does not consider the phenomenon of diffraction “on account of its complicated 
nature”. Plank goes on to say “[i]t will be assumed that the linear dimensions of all parts 
of space considered, as well as the radii of curvature of all surfaces under consideration, 
are large compared with the wavelengths of the rays considered. With this assumption 
we may, without appreciable error, entirely neglect the influence of diffraction caused by 
the bounding surfaces, and everywhere apply the ordinary laws of reflection and 
refraction of light”. Many of the instruments modeled by the Thermal Radiation Group 
do not adhere to this restriction, rather some of the instrument dimensions are on the 
order of the wavelength of the entering radiation. In such cases, general methods used to 
model radiation heat transfer must be supplemented with some means of dealing with 
diffraction in the development of the complete radiative model of an instrument. 

Much of the work done in the Thermal Radiation Group involves the development of 
radiative models of instruments using the Monte-Carlo ray-trace (MCRT) method. This 
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is a statistical approach which simulates the behavior of radiant energy as it enters and 
passes through an instrument. A complete description of this method can be found in 
Walkup [1996], To model radiant energy as it enters an instrument in the MCRT 
environment, the entering energy is divided into a large number of discrete energy 
bundles. These bundles are fired uniformly and randomly into the aperture. For 
example, for the case of a rectangular aperture with comer coordinates of (Xo,Y 0 ) and 
width and length (W,L), illustrated in Figure 4.1, a given energy bundle will enter at 
location X=X<>+R x L, Y=Y 0 +R y W, where R x and R y are random numbers uniformly 
distributed between 0 and 1. Up to this point, radiative models created in the Thermal 
Radiation Group ignored diffraction effects and the entering energy bundles were 
directed parallel into the instrument and treated as collimated energy. It is desirable to 
develop a means of modeling the true behavior of this radiation as it enters and passes 
through the instrument. This proper'modeling of the diffraction of radiant energy as it 
enters radiometric instruments is addressed in this chapter by presenting two candidate 
modeling approaches. 



Figure 4.1 Illustration of physical parameters important to the description of the 

statistical model of diffraction. 
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4.1 Model 1: The statistical approach 

The statistical approach to modeling diffraction was originally published by Heinisch and 
Chou [1971], and later modified and republished by Likeness [1977]. This approach 
treats diffraction of energy as it enters through an aperture in a statistical manner, and is 
designed for application in the Monte-Carlo ray-trace environment. It is based on 
Heisenberg’s uncertainty principle and the particle theory of light. Neither of the two 
references that describe this approach includes many details of its implementation. The 
details necessary to successfully implement this approach, as developed by the current 
author are presented. Note that these details may differ from those used by the original 
developers of this method. 

4.1.1 Appropriate application of the statistical approach 

Rays are directed into the aperture as_ appropriate in the Monte-Carlo environment, and 
the distances 5i and 82 from the point of entry to the aperture edges are calculated for 
each entering ray, as illustrated in Figure 4.1. Experience tells us that entering rays will 
be diffracted; however, the extent to which this diffraction occurs depends upon the 
wavelength of the entering radiation relative to the dimensions of the aperture through 
which it enters. When modeling the diffraction of rays as they pass through an aperture 
which is much larger than the wavelength of the entering radiation, only rays in a limited 
region, or no rays at all, should be diffracted. In some such cases it will be acceptable to 
neglect diffraction altogether and use an ideal ray approximation. Wyatt [1987] suggests 
that this approximation is only acceptable when Aq > 1 0.0 (where A£, is given by equation 
3.1). However, diffraction patterns shown in Figure 3.4 indicate that this rule may be too 
restrictive. Regardless, when modeling a case in which the entering wavelength is much 
shorter than the dimensions of the aperture through which it enters, but where the ideal 
ray approximation would not be appropriate, Likeness [1977] suggests that the statistical 
diffraction model be applied to within 50X to 500X from the edge, while Morbey and 
Hutchings [1993] suggest application to within 100 X. 

If the dimensions of the aperture are on the order of the wavelength of the entering 
energy, all entering rays must be diffracted. The statistical approach treats diffraction as 
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an edge effect, and the question arises as to which edge (or if both edges) of the aperture 
should affect the entering energy bundle. A method of dealing with such cases 
appropriately is needed. Four methods were investigated and are presented in this 
chapter. However from this point forward, the parameter 5 is used as a general parameter 
that takes on both values, 8i and 82, depending on the ray that is being traced and the 
method being used. 


4.1.2 Description of the statistical approach 

This approach is based in part on the Heisenburg uncertainty principle, which states that 
when making simultaneous measurements of a particle’s momentum, p, and position, x, 

the product of the uncertainties in these two measurements cannot be less than — , where 

2 

- *"■ 

h = — , and h is Planck’s constant (6.63 x 10‘ 34 J.s). 

2 n 

That is, 


Ap x 



(4.1) 


where A indicates the amount of uncertainty in a particular measurement. Since a 
particle’s position is described in three dimensions, we can write the uncertainty principle 
in three ways: 


and 


h 

Ap x x Ax >—, 

(4.2) 

h 

Ap y X Ay >-, 

(4.3) 

, h 

Ap z x Az > — . 

2 

(4.4) 


We model Heisenburg ’s uncertainty principle as an equality, 

A . h h 

A P y x Ay = — = — , 

2 4 7t 


(4.5) 


and take the uncertainty in position (in the y direction) of an entering energy bundle to be 
8, the distance from the point of entry to the aperture edge. Thus 
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Ap y = 


4§7c (4-6) 

Referring to the construction of Figure 4.2, the uncertainty in momentum in the y 
direction is interpreted in terms of an angle of diffraction. We define <ji* as the maximum 
angle of diffraction where the angle of diffraction, <J> d , is less than or equal to + <j>*. Thus 
<|>* defines the range of values that <J)d can take on. 



Figure 4.2 Illustration of possible spread of diffraction angles in the application of 

the statistical approach to modeling diffraction. 

Using this construction we can write 


tan(<f>*) = 


Apy 


(4.7) 


li 2n 

Using the equalities pz = — , the wave number, k = — , and equation 4.6, and substituting 

A. A. 

into equation 4.7 yields 


tan(<|>*) = 

P* 


1 

2k5 ’ 


(4.8) 


or 


4>* = tan 1 


(—) 

v 2k5 ; ' 


(4.9) 
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Likeness proposed the use of a normal distribution as the probability density function of 
the diffracted energy, where the most probable path of an entering energy bundle is 

( i 

straight ahead (mean, p=0), and the standard deviation is taken to be <j>* = tan' 1 

v2k5 


Thus 


p(<j)d) 


1 

—= exp 

V2 n <(>* 



(-oo < <j> d < oo) (4.10) 


is the probability that a ray will be diffracted at an angle <j> d from its original line of travel 
due to the influence of a single aperture edge. In order to make use of this probability 
density function in a Monte-Carlo ray-trace environment, we integrate the probability 
density function with respect to <f» d from -oo to <|> d to obtain the cumulative distribution 
function (c.d.f). A c.d.f., P x (x), has the useful property that “if x is a random number, and 
we have the c.d.f, P x (x), and if P x (x) is continuous, the random variable Y produced by 
the transformation Y=P x (x) has a uniform distribution over the interval (0,1)” [Gibbons 
and Chakraborti, 1992]. For a more detailed discussion of this c.d.f. property, see 
Walkup [1996]. We must interpret this integral in terms of the physical limits, where the 
range of angles that <(> d can subtend after entering the aperture is limited to 



Figure 4.3 Illustration of approaching and diffracted ray angles. 
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We must therefore do a substitution such that when |<|>d| > — , p(<f><j) = 0 . Such is the case 


when evaluating p(tan(<J>d)) as illustrated in Figure 4.4. 



<j>, radians 


Figure 4.4 Probability density function of tan(<j>d) 

The c.d.f., which ranges uniformly between (0-1), becomes 


P[tan(4>d)] = j"p(tan<|)'d) d(tan<{»'d) = Random number uniformly distributed (4. 1 1) 




between 0 and 1 , where 


p[tan(4>'d)] = y— ! exp 

V 2n <()* 


tan(<J>'d) 


(4.12) 


Thus 


pMw!= 72^J exp 


^tan(<j)'d)^ 2 
v ♦* j 


d(tan(<j»'d) . 


(4.13) 


-oO 
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Making the substitution 


a 2 = 


2 (< j >*) 2 

there results 


1 1 
, a = 


V2 <t>*’ 


P[tan(<}>d)] = -j= fexp[-a 2 tan 2 (<j)'d)] d[tan(<j>'d)]. 

V27t 


A second substitution, u = a tank’d) , thus du = ad[tan(<|>'d)] yields 

i a(lan^d) , 

P[tan(<J)d)] = -j== — Jexp(- u 2 )— = R* , 

V27T <p* a 

where Rj, is a random number uniformly distributed between 0 and 1 . 
multiplying both sides of equation 4. 1 5 by 2, there results 

a (tan +4) 

2 V 2 <j>*aR* =>-^= jexp(-u 2 ) du. 


Substituting a = 


V2 


■ on the left-hand side yields 


a(tan +d) 


2 = fexp(-u 2 ) du. 

v« -1 

Equation 4. 1 7 can be rewritten 


a(tan fa) 


2R* = —f= fexp(-u 2 ) du + —j= |exp(-u 2 ) du 

V7i.i VTt 0 J 


or 


a (tan fa) 


2R* = — 7 = fexp(-u 2 ) du + —?= fexp(-u 2 ) du. 

Vn 0 0 

Then invoking the definition and properties of the error function, 

2 x 


erf(x) = -= fexp(-ri 2 ) drj , erf(-x) = -erf(x) , erf(oo) = 1 , 
vrc 0 


erf(0) 


we find that 


(4.14) 

(4.15) 
Rearranging and 

(4.16) 

(4.17) 

(4.18) 

(4.19) 

= 0, 


2R> = -erf(-co) + erf[a tan( <t>d)] = l + erf[a tan(<j>d)]. 


(4.20) 
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Finally, substituting a = -j= yields 

V2 <j>* 


2R+ - 1 = erf 


tan(<j>d) 

V2 <j>* 


(4.21) 


Equation 4.21 is used to determine the angle at which the entering ray is diffracted. An 
additional programming step may be required in order to arrive at a value for <t> d - If using 
a programming language without a built-in error function (such as FORTRAN), an 
infinite series approximating the error function and the root-finding bisection method can 
be used to find the value of x when erf(x) is known. This infinite series, given by 



_xl + J cl__x^ + ^ 
3-1! 5-2! 7-3! + J’ 


(4.22) 


when expanded to eleven terms approximates the error function quite well when erf(x) < 
0.992869. Additional terms do not provide much improvement, as shown in Figure 4.5. 
If erf(x) is between 0.992869 and 1, this infinite series deviates from erf(x), and x can be 
approximated as 2 (or any other number greater than 2). The subroutine used is provided 
in Appendix A. Note that the results obtained using this subroutine were compared to 
results found using FORTRAN Powerstation, which includes a built-in error function. 
The resulting diffraction patterns were identical, thus it can be concluded that either the 
necessary approximations had no effect on the results, or that FORTRAN Powerstation 
makes the same approximations. 
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Figure 4.5 Approximating the error function with an infinite series with an increasing 
number of terms. 
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4.1.3 Determining the diffraction pattern 

As previously noted, for cases in which the wavelength of the entering radiation is of the 
same order as the dimensions of the aperture through which rays enter, each entering ray 
must be diffracted in some way. The angle of diffraction depends on the method chosen 
and the point of ray entry. After being diffracted in some way, the energy bundle 
continues until it is intercepted by the observation screen, and this point of interception is 
determined. In order to keep track of the spatial distribution of rays as they arrive at the 
observation screen illustrated in Figure 4.1, the screen is divided into a large number of 
strips in the y direction, serving as discrete “bins”. Each time that a diffracted ray arrives 
at one of these bins, a counter is incremented for that particular bin. This provides the 
distribution of diffracted energy along the y direction. The number of rays arriving in 
each bin is normalized with respect to the bin with the most arriving rays. This 
normalized curve is then scaled such that the area beneath the resulting curve is 
approximately the same as that beneath the curve of the analytical solution (here this has 
been done “by eye”). This was similar to the approach used by Heinsch and Chou 
[1971], as they report using a scaling factor, f, which is determined by trial and error to 
give the best possible match with the analytical solution. The statistical approach does 
not keep track of phase, thus the results obtained do not model the side fringes of the 
analytical interference pattern, rather they form a single smooth curve that averages these 
fringes. Keeping track of phase would not result in a diffraction pattern with the details 
of the side fringes, for reasons that will become clear after the presentation of the second 
model. These results are adequate for most modeling efforts as it will usually suffice to 
know the approximate distribution of diffracted energy, which can be modeled quite well. 
This approach is applicable to both Fraunhofer and Fresnel diffraction. 

4.1.4 Modeling diffraction in a practical example 
4.1.4.1 GERB linear-array cavity detector 

A project that has recently been the topic of much research in the Thermal Radiation 
Group involves a thermal radiation detector originally designed for GERB (Geostationary 
Earth Radiation Budget) [Mahan and Langley, 1996; Weckmann, 1997; Mahan et al.. 
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1998; Sanchez, 1998; Sorensen, 1998]. This detector consists of a linear-array detector 
housed within the wedge-shaped, mirrored cavity shown in Figure 4.6. Radiant energy 
from an Earth scene enters the cavity through a 60 pm wide slit. The detector is used to 
measure the arriving radiation in the wavelength interval between from 0.32 and 40 pm. 
This represents a case in which neglecting diffraction effects in the radiative model may 



The first step in assessing diffraction effects in the cavity detector involves determining 
the regime into which the occurring diffraction falls. We imagine an xy plane passing 
through the thermopile detector along its top edge, and use A5, to determine the 
diffraction regime at this plane. We see that when the entering wavelength (X) is 0.4 pm, 
A£ = 17.3. Likewise, when A. =4.0 pm, A^ = 5.5 and when X=40 pm, A£ = 1.7. Since A4 
>1.0 over the total range of wavelengths to be measured, diffraction within the cavity will 
always occur in the Fresnel region. Note that if we were to study diffraction on an xy 
plane passing through the bottom edge of the detector, Fresnel diffraction would still 
prevail over the measured wavelengths. Figure 4.7 illustrates the changing diffraction 
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pattern which will occur on the linear-array detector as radiant energy of three 
wavelengths enters the GERB cavity detector. 



Figure 4.7 Diffraction patterns expected at the top edge of the linear-array detector 
over the range of wavelengths of interest. 

4.1.4.2 Application of the statistical approach 

The statistical approach was applied to model the diffraction of energy as it enters the 
GERB detector geometry. A radiative model of the GERB cavity detector which ignores 
diffraction effects already exists, so the effort reported here serves as a final step in the 
development of a complete radiative model. 

The diffraction patterns obtained from the application of the statistical approach are 
compared to the closed-form analytical diffraction patterns in order to validate the use of 
the statistical method. This is also a study of an intermediate step within the statistical 
method. Recall that as an energy bundle enters a slit aperture, the distances from the 
point of entry to each edge, 8| and 82, are calculated. It was previously stated that we 
take the unknown in the entering energy bundle’s position as 8. We must establish which 
8 is to be used; that is, we must decide which edge is assumed to diffract the entering 
energy bundle when the slit is on the order of the wavelength of the entering energy. In 
cases such as the aperture of the cavity detector in Figure 4.6, which is essentially an 
infinite slit, it is expected that both edges should affect the momentum of the entering 
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energy bundle in the same direction (in the present case, the y direction). Neither of the 
two references cited address this situation, so four possible approaches were tested. 

4.1.4.3 Four methods of application of the statistical model 

We have developed four methods for using the statistical model to describe diffraction. 
The first method involves summing the angles of diffraction due to both edges of the slit 
as they influence a single ray. As a ray enters, the angle of diffraction due to each edge is 
determined using equation 4.21, one independent of the other, and the two diffraction 
angles are algebraically summed to give the final diffraction angle of the entering ray. 

The second method involves splitting the rays as they enter through the aperture. As a 
single ray enters, this single ray becomes two rays. One ray is diffracted by the first 
edge, located 8 i from the point of entry. The other ray is diffracted by the opposite edge, 
82 from the point of entry. 

The third method involves allowing the side that is closest to the point of entry to diffract 
the entering ray, ignoring the other side. 

The final method involves defining the angle of diffraction of an entering ray as the angle 
at which the ray is diffracted due to the influence of the edge that is furthest from the 
point of entry, ignoring the closest edge. 

All four of these approaches were implemented for five different cases involving the 
GERB cavity detector of Figure 4.6. The angle of approach <|>j (illustrated in Figure 4.3) 
was zero for all cases studied. These cases involved radiation entering at wavelengths of 
0.4, 4.0, 40.0, 100.0, and 160.0 pm; corresponding to values of of 17.3, 5.5, 1.7, 1 . 1 , 
and 0.86, respectively. Note that the cavity detector is not intended to measure 
wavelengths beyond 40.0 pm. Diffraction was studied at 100 and 160 pm in order to test 
the statistical method in both the Fresnel and Fraunhofer regions, as well as in the 
transition region. The FORTRAN code used to conduct this study is provided in 
Appendix A. 
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4. 1.4.4 Results 

Figures 4.8 through 4.12 provide the results of the study described in Section 4. 1.4.3. 
Figure 4.8 involves modeling diffraction in the extremely near-field. This case 
approximates the diffraction situation that occurs in the GERB cavity-detector when used 
to measure radiation at the shortest measured wavelength of 0.32 pm. In this case, any of 
the four methods yields reasonable results. However, results from the summing angles, 
splitting rays, and closest side approaches approximate the analytical solution particularly 
well and are practically identical. According to the previously described criterion by 
Wyatt [1987], this is the only configuration for which the ideal ray approximation would 
be appropriate. Figure 4.9 also involves diffraction in the very near-field. In this case 
both the summing angles and the closest- side approach yield the best results. Figure 4.10 
presents the results from modeling diffraction in the barely near-field, and again the 
summing angles and closest side approaches work the best. This case models the 
diffraction situation that occurs in the GERB cavity-detector when used to measure 
radiation at the longest measured wavelength of 40 pm. Figure 4.1 1 involves modeling 
diffraction in the transition region. In this case the results obtained from the summing 
angles approach are superior to those obtained from application of all of the other three 
cases. Finally, Figure 4.12 involves diffraction in the far-field. In this case, either the 
summing angles or the closest side approach works reasonably well. 
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50.0 0.0 50.0 -50.0 0.0 50.0 

Distance from the central axis, (am Distance from the central axis, (am 


Closest side Furthest side 



-50.0 0.0 50.0 -50.0 0.0 50.0 

Distance from the central axis, jam Distance from the central axis, (am 


Figure 4.8 Comparison between the results using all four methods of implementing the 
statistical approach to diffraction, and the closed-form analytical description of 
the diffraction of radiation entering the cavity detector when the entering 
wavelength is 0.4 jam, the slit width is 60 (am, and the distance to the screen is 
60 jam. 
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40.0 0.0 40.0 40.0 0.0 40.0 

Distance from the central axis, jam Distance from the central axis, (am 

Closest side Furthest side 



40.0 0.0 40.0 40.0 0.0 40.0 

Distance from the central axis, |am Distance from the central axis, (am 

Figure 4.9 Comparison between the results using all four methods of implementing the 
statistical approach to diffraction, and the closed-form analytical description of 
the diffraction of radiation entering the cavity detector when the entering 
wavelength is 4.0 (am, the slit width is 60 jam, and the distance to the screen is 
60 (am. 
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35.0 0.0 35.0 -35.0 0.0 35.0 

Distance from the central axis, pm Distance from the central axis, pm 


Nearest side Furthest side 



-35.0 0.0 35.0 -35.0 0.0 35.0 

Distance from the central axis, pm Distance from the central axis, (am 


Figure 4.10 Comparison between the results using all four methods of implementing the 
statistical approach to diffraction, and the closed-form analytical description of 
the diffraction of radiation entering the cavity detector when the entering 
wavelength is 40 pm, the slit width is 60 pm, and the distance to the screen is 60 
pm. 
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■80.0 0.0 80.0 -80.0 0.0 80.0 
Distance from the central axis, jam Distance from the central axis, (am 


Closest side Furthest side 



-80.0 0.0 80.0 *80.0 0.0 80.0 
Distance from the central axis, (am Distance from the central axis, (am 


Figure 4.11 Comparison between the results using all four methods of implementing the 
statistical approach to diffraction, and the closed-form analytical description of 
the diffraction of radiation entering the cavity detector when the entering 
wavelength is 100.0 jam, the slit width is 60 pm, and the distance to the screen is 
60 pm. 
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Distance from the central axis, pm 



Distance from the central axis, jam 
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Figure 4.12 Comparison between the results using all four methods of implementing the 
statistical approach to diffraction, and the closed-form analytical description of 
the diffraction of radiation entering the cavity detector when the entering 
wavelength is 160.0 pm, the slit width is 60 pm, and the distance to the screen is 
60 pm. 
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4.1.5 Conclusions: Method of choice in application of statistical model 

Results in Figures 4.8 through 4.12 show that the summing angles approach consistently 
yields good results, while the other three methods only yield acceptable results for certain 
cases. The results obtained while using the summing angles approach serve as validation 
of the statistical approach to the numerical modeling of diffraction. 


4.2 Background for understanding diffraction Model 2 

This presentation serves to provide an understanding of the nature of diffraction, and as 
background for understanding the second diffraction model. It borrows heavily from the 
presentation of Fraunhofer diffraction by an infinite slit by Serway [1994], 

The first step in the determination of-a-diffraction pattern resulting from an infinite slit 
involves dividing the slit into a large number (n) of zones of width (Ay) where Ay=a/n, as 



Each zone is to serve as a source of coherent radiation such that each contributes an 
incremental electric field amplitude, AE, at any point on the observation screen. If this 
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approach were modeled in a ray-trace environment, instead of treating the most probable 
path for an entering energy bundle as straight ahead with a spread of paths following a 
normal distribution as in the statistical method, this method would treat each as having an 
equal probability of going in any direction. If a sufficient number of energy bundles were 
directed into the aperture at a given zone, the result would be an even distribution over 
the entire observation screen due to that zone. 


Diffraction occurs because adjacent areas on an aperture all behave as independent 
sources. The energy leaving a given point on the aperture and arriving at a given point P 
on an observation screen will differ in phase from energy arriving to the same point P 
from another point on the aperture. The difference in path length, 5, of rays coming from 
adjacent zones and arriving at the same point on an observation screen can be 
determined. This difference in path length is indicative of the difference in phase, since 
the phase of a ray is proportional to the distance it travels. Suppose that an infinite slit of 
width a is divided into two halves, as pictured in Figure 4. 14. 



Figure 4.14 Illustration of the determination of the difference in path length traveled 

by rays entering from different halves of an aperture. 


The difference in path length between waves entering the two halves is given by 

8 = — sin0 . 

2 


(4.23) 


In general, if a slit aperture is divided into n parts, the difference in path length between 
adjacent areas is given by 
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(4.24) 


Since the width (Ay) of the strips into which the aperture is divided is given by a/n, 

5 = Ay sin 0. (4.25) 


The incremental electric field amplitudes between adjacent zones are out of phase with 
one another by the amount AB, given by 



(4.26) 


where X is the wavelength of the entering radiation, and 5 is the difference in path length 
traveled by energy bundles leaving adjacent areas and arriving at the same point on the 
screen. Using phasor diagrams to determine the intensity at a given point, the total 
electric field is obtained by summing the contributions from all zones at a given point on 
the observation screen. The chord length, Eep, is taken to be the amplitude of the electric 
field at P. Figure 4.15 (b) illustrates the addition of phasors for an infinite slit divided 
into four areas, shown in 4. 1 5 (a). When the number of divisions on the aperture goes to 
infinity, the phasor diagrams become smooth curves, as illustrated in 4.15 (c). 
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(b) (c) 


Figure 4.15 Determination of the amplitude of the electromagnetic field at a point on 
an observation screen placed in front of an infinite slit aperture. (a) 
Illustration of slit division, (b) phasor diagram construction, and (c) 
smooth curve that the phasor diagram becomes as the number of slit 
divisions goes to infinity. 


Taking the arc length to be E 0 , R to be the radius of curvature, and the total phase angle 
from the top to the bottom of the aperture to be [3, simple geometric relations show that 



Ee/2 

R 


Substituting the fact that the arc length, E 0 , is given by Rp yields 


(4.27) 
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E. = EJ 


sin(p/2) 

. (P/2) . 


(4.28) 


Since intensity is proportional to the square of the electric field, we can then write 

~l2 


Ie = In 


sin(p/2) 

(P/2) 


(4.29) 


where I max is the intensity of the central maximum. Because we are discussing the entire 
slit. Ay = a, thus 8 = a sinG. Substituting into equation 4.26, P in this case is given by 

2Tcasin0 

1 ' ( 4 - 3 °) 


P = 


4.3 Model 2: Application of the modified Huygens-Fresnel principle 
4.3.1 Basic description of Model 2 

The second model proposed for use in the modeling of diffraction in the Monte-Carlo 
ray-trace environment is based on concepts presented in Section 4.2. It involves firing 
rays into an aperture, and modeling each point of ray entry as a source itself. The 
original plan for this approach was to model the point of entry of a ray as a source as 
defined by the Huygens-Fresnel principle, whereby the diffracted ray would have equal 
probability of going in any direction and each ray would carry with it an amplitude of 
unity. The distance traveled by the emitted ray before being intercepted by an 
observation screen was to be determined, and the phase assigned to this ray would be 
proportional to its length of travel. After many rays had been traced, the resulting 
intensity distribution was to be determined using a method that will be described shortly. 
The results obtained from the application of this approach did not agree with theory, as 
the secondary peaks were too high relative to the central maxima. A modified approach 
was then implemented in which the point of entry is modeled as a source which emits a 
ray in any forward direction. The amplitude of the optical field at any point on the 
observation screen is taken to be the superposition of all the rays arriving to that point 
(considering their amplitudes and relative phases). This approach is different from the 
initial approach in that all rays do not carry the same amplitude of unity. Instead, the 
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magnitude of each ray is determined by an obliquity factor which is a function of the 
angle of ray diffraction. 


4.3.2 The obliquity factor 

The obliquity factor serves to properly model the variations in amplitude with angle over 
the surface of the secondary wavelets emanating from a source point, modeled incorrectly 
by the Huygens-Fresnel principle [Hecht and Zajac, 1974], Kirchoff proposed that the 
obliquity factor, K, be given by 


K(<|>i, <J>d) s (cos <j>i + cos <(>d) , 


(4.31) 


where <{>j is the angle of incident radiation and <Jm is the angle of diffraction, as illustrated 
in Figure 4.3 [Mayes and Melton, 1994], Later modified by Rayleigh and Sommerfeld, 
and then by Miller, the suggested obliquity factor became 

K(<j>d) = cos(<j>d) . (4.32) 

For the present purposes, use of either of these factors provides comparable results. 
Figure 4.16 shows the weighting to be placed on an diffracted ray as a function of 
diffraction angle when the incident energy approaches normal to the slit, as suggested by 
Huygens, Kirchoff, Rayleigh and Sommerfeld and Miller. 


Figure 4. 1 7 demonstrates how a source point of radiant entry would be modeled if we 
were to apply the Huygens-Fresnel principle alone. When weighted by the 
Rayleigh/Sommerfeld obliquity factor, the entry point is modeled as illustrated in Figure 
4.18. 



Weighting of ray 
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Huygens Kirchoff Rayleigh-Sommerteld/Miller 


1.2 



Figure 4.16 Weighting to be placed on rays when using two different obliquity factors, or 
none at all. 



Z 


Figure 4.17 Model of a point source of rays entering a slit when employing the 
Huygens-Fresnel principle with no obliquity factor. 
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A 



Figure 4.18 Model of a point source of rays entering a slit when employing the 
Huygens-Fresnel principle with the Rayleigh/Sommerfeld obliquity 
factor applied. 


4.3.3 Application of Model 2 

The implementation of Model 2 in the Monte-Carlo ray-trace environment for the 
prediction of the far-field diffraction pattern involves firing energy bundles uniformly 
and randomly across the aperture towards the observation screen which is divided into a 
large number of bins, as in the statistical model. In this case, however, the angle of 
diffraction is chosen so that all angles are equally probable (diffusely scattered) rather 
than following a normal distribution. The length of the path that each energy bundle 
travels and the angle at which it is diffracted are substituted into a modified form of 
equations 4.26 and 4.29. Taking the amplitude of the electric field of a given energy 
bundle to be 


where 


E* - K(<f>i,<j>d) 


Sin(P) 

P 


(4.33) 


P = 2^— and t is the path length. After a large number of rays has been fired (e.g. one 

X 

million), the final intensity in a given bin on the observation screen is the square of the 
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sum of all of these electric fields due to energy bundles arriving at this bin from all parts 
of the aperture, given by 

L. = (£ E.(bin)] ' . < 434 > 

Determining Ibin for all bins on the observation screen and dividing by the maximum 
intensity arriving in any of the bins provides the normalized intensity along the 
observation screen. This technique was tested on two configurations, chosen so that the 
diffraction would be in the Fraunhofer diffraction regime. 

4.3.4 Results from the application of Model 2 

The results from the application of Model 2 are presented as normalized intensity 
compared with the normalized closed- form analytical solution. Figure 4.19 presents the 
results from application of this approach to an infinite, 0.3-mm-wide slit with energy 
entering at the wavelength of 0.58 pm, and an aperture-to-screen distance of 2 m. 


• Ray trace Analytical 



Figure 4.19 Comparison of results from application of modified Huygens-Fresnel 
principle and analytical solution for far-field diffraction from an infinite slit. 




Normalized intensity 
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Note that this approach produces a pattern with high-frequency oscillations in intensity, 
whose peaks lie within an envelope that matches the closed-form analytical solution quite 
well. 

It is interesting to study the application of Model 2 to other aperture shapes, such as a 
circular aperture. A circular observation screen divided into bins consisting of equal-area 
rings is placed behind this aperture, and the same procedure is followed to obtain a plot 
of intensity with angle from the central axis. The case studied involved a circular 
aperture 0.2 mm in diameter with entering energy of wavelength 9 pm, and aperture to 
observation screen distance of 16 m. 


• Ray trace Analytical 



0.00 0.02 0.04 0.06 0.08 0.10 

Angle a from central axis (radians) 


Figure 4.20 Comparison of results from application of modified Huygens-Fresnel 
principle and analytical solution for far-field diffraction from a circular 
aperture. 
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The resulting intensity pattern shown in Figure 4.20, when compared to the analytical 
solution, is more difficult to interpret than that for the infinite slit. For this case, two 
oscillating patterns, or orders, occur which together form an envelope that approximately 
matches the analytical solution. The oscillating intensity patterns seen in both cases 
occur for reasons that are not clear. For the case of the slit, the resulting pattern is similar 
to that formed by radiation passing through multiple closely-spaced slits [Hecht and 
Zajac, 1974], The agreement between the diffraction predicted by this model, and the 
closed-form analytical solution will occur in certain cases when the Fraunhofer condition 
applies, described in Sections 4.3.5 through 4.3.7. Otherwise, the simplifying 
assumptions break down, and we can no longer expect good results from this application 
of the modified Huygens-Fresnel principle [Hecht and Zajac, 1974], The FORTRAN 
code used to generate these results is provided in Appendix B. 

4.3.5 Limitations of the Huygens-Fresnel principle and of Model 2 

Because Model 2 is based upon the modified Huygens-Fresnel principle, it is important to 
understand the limitations of this principle. As stated by Hecht and Zajac [1974], in 
cases in which the aperture is very large, and the point of observation is far away, the 
Huygens-Fresnel principle should, and does, work very well. However, for cases 
involving a very small aperture, or when the point of observation is in the vicinity of the 
aperture, deviation from the behavior predicted by the Huygens-Fresnel principle should 
be appreciable. Here, the size of the aperture refers to the aperture dimensions relative to 
the wavelength of the entering radiation. These limitations imply that Model 2 may 
properly model diffraction only for configurations where A£ < 1.0, and where the 
aperture width-to-wavelength ratio (a/7.) exceeds some minimum value. In order to 
investigate these limitations, several cases were studied in which the value of A^ and the 
aperture width were held constant, while the ratio of a/X was varied. These cases were 
studied by comparing the closed-form analytical solution of the diffraction pattern with 
the results predicted by Model 2. 

The case study involved modeling diffraction by an infinite, 60 pm-wide slit. The 
wavelength of the entering radiation and the aperture-to-observation screen distance were 
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varied so that was a constant 0.4. Five cases were studied, whereby aJX was assigned 
the values of 4.0, 1.5, 1.0, 0.5, and 0.25. The results from this study are provided in 
Section 4.3.6. 

4.3.6 Results from case study 

Model 2-predicted results shown in Figure 4.21 agree well with the analytical solution 
(the envelope containing the peaks of the oscillating pattern predicted by Model 2 
approximately match the closed-form analytical solution). For this case, the entering 
wavelength is four times smaller than the aperture width and the observation points are 
far from the aperture; as expected, the Huygens-Fresnel principle works well. Figure 
4.22 demonstrates the declining performance of the Huygens-Fresnel principle as the 
wavelength of the entering radiation approaches the width of the slit. Here the central 
maxima is outlined by the peaks of the oscillating pattern predicted by Model 2, but the 
secondary fringes are lost entirely. The results shown in Figure 4.23 are for the case in 
which the entering radiation is of the wavelength equal to the slit width. These results are 
beginning to lose the ability to predict even the shape of the central maxima. As the 
entering wavelength is made larger than the slit width, Model 2 produces highly 
erroneous results, as illustrated in Figures 4.24 and 4.25. 

• Ray trace Analytical 



0 500 1000 1500 

Figure 4.21 Comparison of results from application of modified Huygens-Fresnel 
principle and analytical solution for far-field diffraction from an infinite slit 
aperture for which A£ = 0.4 and a JX = 4.0. 
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• Ray trace 


Analytical 



0 500 1000 1500 


Figure 4.22 Comparison of results from application of modified Huygens-Fresnel 
principle and analytical solution for far-field diffraction from an infinite slit 
aperture for which A£ = 0.4 and a IX = 1 .5. 



0 500 1000 1500 


Figure 4.23 Comparison of results from application of modified Huygens-Fresnel 
principle and analytical solution for far-field diffraction from an infinite slit 
aperture for which = 0.4 and a/X = 1.0 
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• Ray trace Analytical 



o 500 • ' 1000 1500 


Figure 4.24 Comparison of results from application of modified Huygens-Fresnel 
principle and analytical solution for far-field diffraction from an infinite slit 
aperture for which A£, = 0.4 and a /X = 0.5. 



Figure 4.25 Comparison of results from application of modified Huygens-Fresnel 
principle and analytical solution for far-field diffraction from an infinite slit 
aperture for which A£ = 0.4 and a /X = 0.25. 
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4.3.7 Conclusions: Model 2 

The results from the case study described in Section 4.3.6 suggest that Model 2 will 
approximate the diffraction pattern of radiant energy entering an aperture, including the 
secondary maxima only if Ad; < 1 .0, and a A. » 1 .0. Results indicate that as long as a JX > 
1 .0, application of Model 2 will not lead to highly erroneous results, but the details of the 
secondary fringes may be lost. However, if a/ A, < 1.0, Model 2 should not be applied as 
its underlying principles are no longer sound. 

4.4 Conclusions: Model 1 versus Model 2 

The statistical model is the most general choice for the modeling of diffraction in a 
Monte-Carlo environment, as it can approximate interference patterns caused by the 
diffraction of energy for both near and far-field conditions. This model cannot be 
modified to predict the fringes about the- central maxima by keeping track of phase angle 
(P) as is done in the second model. This impossibility is due to the requirement that all 
points on the aperture contribute equally to all bins on the observation screen for the 
formation of an interference pattern. In other words, the number of arriving rays (coming 
from any randomly located point on the aperture) must be the same for each bin, but the 
resulting intensity is given by the sum of the intensity from each individual ray, weighted 
by its phase angle. In the statistical approach, this criterion is not met, as the distribution 
of energy from a given point on the aperture is not diffusely distributed, but follows a 
normal distribution. Model 2 can only be used with good results for a restricted set of 
conditions, otherwise the underlying assumptions break down and it behaves very poorly. 

4.5 Potential future investigations of diffraction models 

The model based upon the modified Huygens-Fresnel principle (Model 2), and the cause 
of the oscillating intensity pattern it predicts, should be further investigated. Other case 
studies similar to the one presented in Section 4.3.6 could be conducted to determine the 
generality of the conclusions in Section 4.3.7. Another unexplored aspect of the 
diffraction of radiant energy by an aperture involves the behavior of radiation as it 
approaches an aperture (prior to entry through the aperture), and the possibility of 
backwards diffraction. The geometric theory of diffraction described in Chapter 3.0 does 
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model both forward and backward diffraction of energy as it approaches an aperture, and 
may be capable of properly modeling this behavior. Useful future work could also 
involve the development of an approach to scale ray-trace results to match the analytical 
curve for the statistical method, as the current method involves simply scaling the results 
“by eye” until the area under the ray-trace curve appears to match that of the analytical 


curve. 



5.0 NEXT-GENERATION INSTRUMENT CONCEPTS 

Studies of the potential next-generation earth radiation budget instrument, 
PERSEPHONE, as described in Chapter 2.0, require the use of a radiative model of the 
CERES and/or modified CERES telescope. The first task of the current effort was to 
develop such a model which is flexible enough to perform the required studies. 
Repeating the research objectives of this thesis topic, the first goal is to determine the 
maximum number of detectors that can be placed in the current CERES telescope without 
the loss of Optical Point Spread Function (OPSF) quality. The other principal objective 
is to investigate the possibility of using hyperbolic mirrors in place of the current 
spherical mirrors in order to maximize the number of channels that can be fitted into the 
current CERES envelope, thus maximizing the science data return. Descriptions of these 
studies, and their results are presented in this chapter. 

5.1 Development of the radiative model of the CERES optics 

Members of the Thermal Radiation Group have developed several radiative models of the 
CERES and ERBE telescopes. Meekins [1990] completed a Monte-Carlo-based 
numerical model to study the optical and radiative characteristics of the ERBE scanning 
radiometer. Bongiovi [1993] later modified Meekins’ code, adding baffles in order to 
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model the CERES instrument. Bongiovi’s code consists of 16,632 lines of commented 
FORTRAN code, and is specific to the CERES geometry. Modification of this code to 
conduct the study at hand would be a formidable task. Instead, a new radiative model 
was developed using a tool which is being developed by Felix Nevarez, a doctoral 
student in the Thermal Radiation Group. This tool consists of a C ++ class library which 
can represent quadric surfaces, such as cylinders, spheres, planes, and cones with 
specified surface absorptivity and specularity ratio. The radiative model of an instrument 
can be developed by modeling each interior instrument surface using one of the library- 
defined surfaces with appropriate dimensions, and stacking these surfaces to model the 
entire instrument. After the complete geometry is described, a Monte-Carlo engine 
module is called, and a user-specified number of energy bundles is directed through the 
instrument. The paths of these energy bundles are determined by following the rules of 
the Monte-Carlo ray-trace method as. it applies to radiation heat transfer. The output 
from the execution of this C code is a file of the distribution of energy bundles arriving 
on a surface of interest. 

5.1.1 Modeling the CERES geometry 

In the present instrument study, the surface of interest is the plane which would contain 
an array of thermistor bolometer detectors, located beneath the precision aperture 
pictured in Figure 5.1 (b). Figure 5.1 (a) demonstrates how the CERES telescope is 
broken down into 46 basic surfaces. The numbering scheme corresponds to the 
numbering used in the code provided in Appendix C. Note that with the use of this new 
tool, the author needed only to write 616 lines of commented C code which calls the C~ 
library functions, and 224 lines of FORTRAN code for post-processing the output files, 
for a total of 840 lines of code. When compared to Bongiovi’s 16,632 lines of 
FORTRAN required to achieve the same task, this new code proves to be a remarkable 
tool in radiometric instrument design. 
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In addition to the ease of performing radiometric analyses, this new ray-trace 
environment is equipped with a useful graphical user interface. After entering the 
telescope geometry, it can be viewed by specifying all of the surfaces as diffuse and 
firing rays into the telescope. Figure 5.2 shows the CERES telescope with the primary 
mirror removed, as produced by this graphical user interface. 



Figure 5.2 CERES telescope geometry, as produced by graphical user interface of the 
new ray-trace environment. 
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5.1.2 Addition of shims for blur circle minimization 

As described in Chapter 3.0, the presence of spherical aberration in an optical system 
causes the spot size at the point of “best focus” to be finite. This spot of finite diameter is 
called the blur circle. It is desired that the location of the minimum blur circle occur at 
the plane containing the precision aperture so that the edges of the Optical Point Spread 
Function demonstrate the steepest drop-off possible (i.e. approach a “top-hat” response). 
The location of the “best focus” is very sensitive to the distance between the primary and 
secondary mirrors. Upon placing the parts between these mirrors, tolerance stack-up 
occurs so that the location of the minimum blur circle will likely not occur at the plane 
containing the precision aperture. The spacing between the mirrors is thus increased by 
the addition of very thin shims until the blur circle is a minimum at the precision 
aperture. This addition of shims was simulated using the computer model of the CERES 
telescope. The distance between the mirrors was gradually increased, and a scattergram 
representing the blur circle was plotted until the distance providing the minimum blur 
circle was found. The change in the image at the focal plane due to slight changes in the 
primary-to-secondary mirror spacing is illustrated in Figure 5.3. The quantity 5 refers to 
the total thickness of the shims added. When 8 = 0.28 mm, the minimum blur circle is 
obtained, as illustrated in Figure 5.3 (b). Note that the structure which supports the 
secondary mirror, the “spider”, is imaged in the defocused images shown in Figures 5.3 
(a) and (c). 


5 = 0.0 


5 = 0.28 mm 


8 = 0.56 mm 



Figure 5.3 Illustration of change in the images at the precision aperture for different 
shim thicknesses, where 8 is the shift due to the addition of shims 
(dimensions in mm), (a) Defocused image, (b) image at best focus where 
the blur circle diameter is a minimum, and (c) defocused image. 
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5.2 Capabilities of the current CERES instrument 

5.2.1 Determination of the Optical Point Spread Function of the current CERES 
instrument 

The OPSF (Optical Point Spread Function) of an instrument shows how collimated 
energy entering from various angles is transmitted through the instrument. The OPSF has 
previously been determined for the CERES telescope with only a single aperture using 
Bongiovi’s model [Priestley, 1997], This OPSF exhibits attenuation at the edges of the 
field due to the finite blur circle. If the blur circle were infinitely small, as in an ideal 
optical system, this edge attenuation would be absent from the OPSF. This ideal OPSF 
would be perfectly flat across the aperture area, and would drop immediately to zero 
beyond the edges. Instead, attenuation occurs before the physical edge of the field stop is 
reached, and radiation arrives at the detector beyond the projection of the area of the field 
stop on the detector. A complete description of the significance of the blur circle is 
presented by Priestley [1997]. 

Equivalence is assumed in the determination of the OPSF of the CERES telescope. This 
assumption implies that regardless of the point at which a given amount of energy arrives 
at a thermistor bolometer detector, the response of the detector is the same. 

In order to determine the OPSF over more than one aperture using the new radiative 
model, collimated radiation was allowed to arrive at angles 0, b (see Figure 5.1 (a)) and 
traced through the instrument. The output file for a given combination of angles consists 
of the x, y coordinates of all energy bundles arriving at the plane containing the precision 
aperture. In order to determine how much of this arriving energy enters one of the 
precision apertures, each output file must be opened and read, line by line. If the 
coordinates of a line of output fall within the area defined by of one of the precision 
apertures, then an energy bundle will reach the detector below, and a counter for that 
aperture is incremented. After all lines of all output files have been read, the results are 
the number of energy bundles arriving to each detector for each set of input angles 0, <J>. 
This is the information needed to construct the OPSF of the instrument. The FORTRAN 
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code used to perform this post processing is provided in Appendix C. Figures 5.4 and 5.5 
show the resulting OPSF of the current CERES instrument with five (5) apertures and 
two (2) apertures, respectively. Section 5.2.3 provides an interpretation of these figures. 
130-0.1 ■0.1 -0.2 □0.2-0.3 aQ.3-0.4 ■0.4-0, 5 ■0.5-0.6 ■0.6-0.7 □ 0-7-0.8 ■0.8-0.9 ■Q.»1 \ 


Scan angle 




Scan angle 

Figure 5.4 OPSF for the current CERES telescope (minus the primary mirror insert) 
with five (5) precision apertures. 
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Figure 5.5 OPSF for the current CERES telescope with two (2) precision apertures. 
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5.2.2 Validation of results 

The minimum blur circle found in the current study is 0.14 mm in diameter, comparing 
relatively well with the 0. 122 mm diameter blur circle reported by TRW [Carman, 1993], 
There is some discrepancy between the OPSF previously determined [Priestley, 1997], 
and the results of the current research effort in that the cutoff at the edges of the previous 
OPSF is not as steep as that of the results presented in this thesis. It is believed that this 
discrepancy occurred because the code used in the previous study did not model the 
telescope at best focus. Figure 5.6 compares results from the current and previous studies 
to those predicted by linear optics, and experimental data borrowed from Priestley 
[1997], The curve predicted by linear optics was constructed by moving a blur circle of 
0.14 mm in diameter across the detector and plotting the fraction of the blur circle area 
which would fall within the precision aperture. Figure 5.7 shows very good agreement 
between the results predicted by linear optics and current results, and reasonable 
agreement with experimental data. 
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Figure 5.6 Comparison of current results with previous results, linear optics, and 
experimental data. 
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5.2,3 Conclusions: Use of existing CERES instrument with spherical mirrors 

Figure 5.4 illustrates that the placement of five detectors in the current CERES 
instrument does not yield acceptable performance. A flat response across all detectors is 
desired, but a drop off in the response occurs at the two detectors furthest from the central 
axis. Note that the only modification to the current CERES telescope involved the 
removal of the primary mirror insert. Figure 5.5 shows that the placement of two 
detectors within the current CERES telescope does yield acceptable results, as the 
response is flat over the two detectors. It can be concluded that the capability of the 
current CERES instrument can be doubled by the addition of another detector in each 
telescope without sacrificing the quality of the OPSF. As previously mentioned, this 
modification would pose other challenges such as the possibility of optical cross-talk 
between channels. 

5.3 Replacement of spherical mirrors with hyperbolic mirrors 
5.3.1 Optical prescription for hyperbolic mirrors 

In initial discussions of the potential next-generation design concept, PERSEPHONE, 
NASA engineers hypothesized that the replacement of spherical mirrors with hyperbolic 
mirrors within the same telescope could yield an instrument with good throughput over a 
significantly larger field of view, and with the same or smaller size blur circle. The basic 
parameters and layout of the current and potential future mirror systems are illustrated in 
Figure 5.7. The primary-to-secondary mirror spacing (vertex-to-vertex) is to be held the 
same, and only the curvatures and conic constants of the mirrors are to be changed. 
NASA engineers performed a preliminary study to determine the set of hyperbolic 
mirrors which met the required spacing and restricted mirror diameters that would yield 
the best performance. The constraints were entered into a commercial ray-trace code 
which iterates until finding the optimal combination resulting in the minimum blur circle. 
The recommended prescription is provided in Table 5.1. 
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Spherical secondary Spherical primary 
mirror mirror 



(a) 


Hyperbolic Hyperbolic primary 



Figure 5.7 (a) Illustration of current CERES spherical mirrors (Modified Cassegrain), 

and (b) illustration of the new optical prescription with hyperbolic mirrors 
(Ritchey Cretian Cassegrain) (dimensions in mm). 


Table 5.1 Optical prescription for hyperbolic mirrors. 



Primary Mirror 

Secondary Mirror 

r v (vertex radius) 

36.042 mm 

32.284 mm 

Max R (maximum mirror radius) 

10.0 mm 

5.0 mm 

C.C. (conic constant) 

-1.3329 

-22.4729 

Primary-to-secondary (vertex-to-vertex) spacing: 

10.89 mm 

Primary vertex-to-detector plane spacing: 2.0 mm 









Katherine L. Coffey 


Chapter 5.0 Next generation instrument concepts 


76 


This combination of mirrors is said to yield blur circle sizes of 0.058 mm diameter on- 
axis, and 0.080 mm at 2.34 mm off-axis. The first objective in conducting this study was 
to use the new ray-trace tool to duplicate these results using only the mirrors. Upon 
successfully modeling this mirror combination, the next objective was to determine the 
Optical Point Spread Function of the hyperbolic mirrors and compare results with that 
obtained for the spherical mirrors. 


In order to model the desired combination of mirrors, the parameters in Table 5.1 had to 
be converted into input required by the new MCRT environment. Details of this 
conversion are provided for two reasons: (1) in order to document the approach for 
future researchers, and (2) in order to clarify the appropriate use of several equations 
found in [Walkup, 1993], a reference commonly used by members of the Thermal 
Radiation Group. - - 


5.3.2 Conversion of known parameters to required parameters 

In order to specify hyperbolic mirrors using the new MCRT environment, the user must 
supply the parameters a, b, and c h , parameters in the standard equation for a hyperboloid 
given by 


f ,, V 


z 

Vch; 


y 

Vb ) 


" I - 1 - 0 ' 


(5.1) 


where a-b for symmetrical optics. The parameters a and ch are illustrated in Figure 5.8. 
Note that a is not the outer radius of the mirror slice, as stated by Walkup [1993], 



Figure 5.8 Illustration of parameters required for entry of hyperbolic mirrors into new 
MCRT environment. 
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Using the two expressions relating the parameters for a hyperboloid. 


and 



C.C. = - 


„2 . „ 2 
a + Ch 


Ch 


(5.2) 


(5.3) 


we have two equations and two unknowns. Therefore the parameters Ch and a can be 
determined. Because these are symmetrical optics, a = b. Finally, the dimension Zc P , is a 
required parameter, and could be found by substituting b = 0, a = Max R, and the value 
of Ch previously determined into equation 5.1. Using these parameters, the prescribed 
hyperbolic mirror combination was modeled with the new MCRT environment. 


5.3.3 Results 

The only results provided from the preliminary study conducted by NASA involved the 
size of the blur circle at several points on the focal plane. The central blur circle diameter 
was stated to be 0.058 mm when the vertex-to-vertex distance between mirrors, Az, was 
10.89. A similar result was obtained when the prescribed mirror combination was 
modeled using the new ray-trace tool, but the shape of this blur circle indicated that the 
mirrors are not at their best focus when Az = 10.89 mm. Increasing Az slightly resulted 
in a much smaller blur circle, as shown in Figure 5.9. The configuration yielding the 
smallest blur circle, where Az = 10.93 mm, was used for the remainder of the studies 
conducted. 
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A z = 10.89 mm 


A z = 10.93 mm 


A z = 10.96 mm 



Figure 5.9 Illustration of blur circle minimization (where Az indicates the distance 
between the primary and secondary mirrors) (dimensions in mm), (a) 
Defocused image, (b) image at best focus where the blur circle diameter is 
a minimum, and (c) defocused image. 


The Optical Point Spread Function (OPSF) for the hyperbolic mirror combination was 
produced using the new MCRT environment. These results are shown in Figure 5.10 (b) 
compared to an OPSF generated for the spherical mirrors used in CERES, removed from 
the telescope (Figure 5.10 (a)). Only one half of the OPSFs are shown, as symmetry 
about the central axis is expected. Results presented in Figures 5.10 (a) and (b) serve to 
compare the performance of the spherical and hyperbolic mirror combinations alone, 
independent of the influence of the telescope geometry. Note that these OPSFs differ 
from those that would be obtained if the mirrors were placed within the telescope. 

It is important to realize that the outer diameters of the hyperbolic mirrors are larger than 
those of the spherical mirrors; thus, some slight changes in the telescope would be 
required in order to use the recommended hyperbolic mirrors. NASA stated that this 
increase in diameter was a requirement for proper performance of the hyperbolic mirrors, 
and hyperbolic mirrors of the same diameters as the current spherical mirrors would not 
provide acceptable results. 
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Figure 5.10 (a) 0PSF of the current CERES spherical mirrors only, and 

(b) 0PSF of the prescribed hyperbolic mirrors only (profile view). 
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Figure 5.11 (a) OPSF of the current CERES spherical mirrors only, and 

(b) OPSF of the prescribed hyperbolic mirrors only (relief view). 
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5.3.4 Conclusions: Hyperbolic mirrors 

Figures 5.10 and 5.1 1 indicate that the hyperbolic mirrors provide the ability to produce a 
uniform flux over a larger area at the focal plane. Comparison of the OPSF of the 
spherical mirrors alone, shown in Figures 5.10 (a) and 5.11 (a), with that of the CERES 
telescope with the spherical mirrors (Figure 5.4) indicates that the drop-off in response is 
largely due to the manner in which radiation is transmitted through the telescope. (Here 
we are careful to distinguish between the full CERES telescope model; which includes 
baffles, secondary-mirror support struts, and other structural members; and an ideal 
optical system having the same prescription but consisting of only the two mirrors and 
the precision apertures). The drop-off due to the real telescope geometry would be 
expected in an instrument containing the hyperbolic mirrors as well, thus results do not 
show that a change in the mirrors will permit the use of five detectors in a given 
telescope. What is demonstrated is the improved performance of hyperbolic mirrors over 
spherical mirrors. Note, however that the current spherical mirror combination involves 
45 percent obscuration while the hyperbolic mirror combination involves 55 percent 
obscuration. The result of such an increase in obscuration is less energy throughput and 
more problems with diffraction. 

5.4 Potential future investigations 

Future studies could involve the placement of the hyperbolic mirrors into an 
appropriately modified CERES telescope, and the determination of the resulting Optical 
Point Spread Function. The new MCRT environment could be used to study the effect of 
slight modifications of the interior telescope surfaces on the throughput to the detectors. 

Another potential study could involve a different detector arrangement within the 
telescope. For example, instead of inserting extra detectors in a row, the detectors could 
be placed in some alternative arrangement which provides an optimal throughput to all 
detectors. One potential rearrangement for the insertion of four precision apertures is 
illustrated in Figure 5.12. 
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Figure 5.12 (a) Arrangement of precision apertures that has been considered, and 

(b) a potential arrangement of precision apertures for future study. 



6.0 CONCLUSIONS AND RECOMMENDATIONS 

6.1 Conclusions: modeling diffraction in the MCRT environment 

Two diffraction models for use in the Monte-Carlo ray-trace (MCRT) environment have 
been described and tested: Model 1, a statistical approach, and Model 2, application of 
the modified Huygens-Fresnel principle. The derivation and application of these models, 
and the results that they predict are presented in detail in Chapter 4.0. 

The concepts upon which Model 1, the statistical approach, is based were originally 
proposed by Heinisch and Chou in 1971 . The details of the implementation of this model 
have evidently never been thoroughly documented in the public domain, so they were 
rediscovered and documented in this thesis. The statistical approach is useful for 
predicting diffraction in both the Fraunhofer and Fresnel regimes, but does not keep track 
of phase and thus cannot exactly predict the details, such as the secondary maxima of the 
expected diffraction patterns. Intensity distributions predicted by Model 1 are smooth 
curves, which approximate the diffraction patterns predicted by theory. The statistical 
approach was tested by applying it to a practical example involving the diffraction of 
radiant energy as it passes through the infinite-slit aperture of a cavity detector developed 
within the Thermal Radiation Group. This case involves the diffraction of energy along 
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the same coordinate direction due to the two edges of the slit, a situation not discussed in 
the cited references. A study was conducted to determine the best approach to model this 
type of situation and it was determined that the final angle of diffraction of an entering 
energy bundle is best described by the algebraic sum of the angles of diffraction due to 
each edge. 

Model 2 involves the application of the Huygens-Fresnel principle, modified by a 
correcting obliquity factor. This model is capable of predicting the diffraction pattern, 
including the secondary maxima, for a limited range of Fraunhofer diffraction 
configurations. A case study involving diffraction by an infinite slit was conducted to 
determine the general range of applicability of this model. The results from this study 
suggest that Model 2 will approximate the diffraction pattern of radiant energy entering 
an aperture, including the secondary .maxima, only if the ratio of the slit width to the 
wavelength of the entering energy, a /X, is much less than one. Results indicate that as 
long as a /X exceeds one, the application of Model 2 will not lead to highly erroneous 
results but the details of the secondary fringes may be lost. However, if a/X is less than 
one. Model 2 should not be applied as its underlying principles are no longer sound. 

6.2 Conclusions: CERES follow-on instrument 

A radiative model of the current CERES telescope was developed using a new Monte- 
Carlo ray-trace environment being developed by a doctoral student in the Thermal 
Radiation Group. This model was used to study the feasibility of partitioning the current 
CERES telescope so that it serves multiple detectors. Also studied was the replacement 
of the spherical mirrors currently used in CERES with hyperbolic mirrors in order to 
achieve acceptable radiative throughput over a larger field of view. 

The Optical Point Spread Function (OPSF) at a single detector placed on the CERES 
optical axis was determined using this new radiative model. This result was used to 
benchmark the new ray-trace environment and to validate the new code describing the 
CERES instrument by comparing it to the OPSF previously obtained by members of the 
TRG. 
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Results from this study indicate that the radiation throughput to two detectors placed in 
the CERES telescope is similar to that arriving at the single detector that lies on the 
optical axis in the current CERES instrument. However, if more than two detectors are 
placed into the telescope, the throughput drops off, resulting in an unacceptable OPSF at 
all but the central detector. It was also determined that hyperbolic mirrors do achieve 
acceptable radiative throughput over a larger field of view than do spherical mirrors, and 
their use is suggested if more than two detectors are to be placed in each telescope. 

These results are presented in Chapter 5.0. 

6.3 Potential future investigations of diffraction models 

Useful future work could involve the development of an approach to scale ray-trace 
results to match the analytical curve for the statistical method, as the current method 
involves simply scaling the results “by eye” until the area under the ray-trace curve 
appears to match that of the analytical curve. The model based upon the modified 
Huygens-Fresnel principle (Model 2), and the cause of the oscillating intensity pattern it 
predicts, should be further investigated. Other case studies similar to the one presented in 
Section 4.3.6 could be conducted to determine the generality of the restrictions placed on 
the applicability of this model. Another unexplored aspect of the diffraction of radiant 
energy by an aperture involves the behavior of radiation as it approaches an aperture 
(prior to entry through the aperture), and the possibility of backwards diffraction. The 
geometric theory of diffraction described in Chapter 3.0 does model both forward and 
backward diffraction of energy as it approaches an aperture, and may be capable of 
properly modeling this behavior. 

6.4 Potential future investigations of CERES follow-on instrument 

Future studies could involve the placement of the hyperbolic mirrors into an 
appropriately modified CERES telescope, and the determination of the resulting Optical 
Point Spread Function. The new MCRT environment could be used to study the effect of 
slight modifications of the interior telescope surfaces on the throughput to the detectors. 
Another potential study could involve a different detector arrangement within the 
telescope. For example, instead of inserting extra detectors in a row, the detectors could 
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be placed in some alternative arrangement, which provides an optimal throughput to all 
detectors. 
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This program , written in FORTRAN, demonstrates the application of the statistical approach to the 
modeling of diffraction. It predicts the spatial distribution of diffracted energy as it passes through a single 
infinite slit, and arrives at an observation screen some distance away. The following describes the 
variables used throughout the program. This particular version is for the splitting rays approach. 

* NY is an array that divides the observation screen into strips. Each element in this array serves as a 

* '‘bin' 1 which keeps track of the number of bundles reaching that strip. 

* LAMDA is the wavelength of the entering monochromatic radiation. 

* W is the width of the rectangular slit. 

* Z is the distance from the aperture to the observation screen. 

* R is the change in vertical location that the entering energy bundle undergoes between its point of entry 

* and its arrival to the observation screen. 

* DIST is the total distance traveled by energy bundle from aperture before being intercepted by the 

* observation screen. 

* YO is the coordinate of the bottom edge of the rectangular aperture. 

* Y is the coordinate of the entering energy bundle in the plane containing the aperture. 

* YSCREEN is the coordinate of the energy bundle when it strikes the observation plane. Note that the 

* YSCREEN coordinate is such that YSCREEN=0 occurs at the aperture center. 

* YSCMIN, YSCMAX determine the minimum and maximum values of Y for which the number of 

* energy bundles striking at the observation screen will be recorded. 

* H is the number of strips into which the observation screen will be divided. 

* FNCREM is the width of each of the strips on the observation screen. 

* NUMRAYS is the number of energy bundles to be directed from the aperture to the observation screen. 

* DEL 1 ,DEL2 are the distances from the point of entry to the two aperture edges. 

* RANMAR calls a random number from the subroutine rmarin(ij, kl). 

* XVAL(ERFX) is the subroutine which determines the values of x when erf(x) is known, and returns this 

* value to the main program. 

* ERFX1,ERFX2 are the arguments sent to XV AL for the determination of x where erf(x) is known. 

* SDl,SD2 are the standard deviations for the diffraction angles calculated for each of the two edges. 

* K is the wavenumber, given by 2*PI/LAMDA. 

* PHI1,PHI2 are the two diffraction angles due to each aperture side. 

* YDIV is the variable used to determine the "bin" into which to store energy bundles incident to the 

* observation screen. 
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PROGRAM DIFFRACT1 

CS NOEXTENSIONS NOWARNINGS 

* Initialize arrays, Note that array NY is of size H, where H is the number of strips into which 

* observation screen is divided. 

DIMENSION NY(1000) 

* Initialize variables used in program 

DOUBLE PRECISION LAMDA, W,Z, Y,R1 , Y 0,DEL 1 ,DEL2,RANMAR,ARG2 
•PI.ERFX 1 ,ERFX2,SD 1 ,SD2,K,PHI 1 ,PHI2,PHIDIFF 1 .YSCREEN 1 ,ARG1, 
*YSCMIN,YSCMAX,INCREM,NY,DISTl,PHIDIFF2,R2,DIST2,YSCREEN2 

REAL XVAL 

INTEGER N,YDIV,H,NUMRAYS 

DATA LAMDA, W,Z,PI/1 00.0,60.0,60.25,3. 14 1 593/ 

H=1000 

YO=-W/2 

YSCMAX=80.0 

YSCMIN=-80.0 

INCREM=(YSCMAX-YSCMIN)/H 

NUMRAYS=2000000 

K={2 * PI)/L AMD A 

* Open file that will store output from execution of code. 

OPEN( 1 5,FILE=’sr 1 00',STATUS='OLD') 

Initialize the values in the matrix containing the number of energy bundles arriving at screen, which is 

* divided into H strips. 

DOS 1=1, H 
NY(I)=0 
5 CONTINUE 

* Assign seeds for random number generator. 

1=1 

J=3 

CALL RMARIN(I.J) 

Begin Monte-Carlo solution. This loop causes NUMRAYS rays to be fired in randomly and uniformly 

* from aperture. 

DO 10 J=l, NUMRAYS 

* Calculate point of entry of the current energy bundle. 

Y=YO + RANMAR()*W 

* Calculate distance from point of entry of the current energy bundle and the two aperture edges. 
DELl=(W/2)+ABS(Y) 
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DEL2=(W/2)-ABS(Y) 


* Calculate the standard deviation for distribution of diffraction angles. 

SD 1 =ATAN( 1 ,0/(2.0*DEL 1 * K)) 

SD2=ATAN( 1 ,0/(2.0* DEL2* K)) 

* Calculate the angle of diffraction caused by each of the aperture edges. The subprogram XVAL must 

* be called in order to determine the value of x, where erf(x) is known. 

ERfXl= 2.0*RANMAR()-1.0 
ARG1=XVAL(ERFX1) 

PHI1=ATAN(ARG1*SQRT(2.0)*SD1) 

* Calculate the angle of diffraction caused by the other aperture edge. 

ERFX2= 2.0*RANMAR()-1.0 
ARG2=XVAL(ERFX2) 

PHI2=ATAN(ARG2*SQRT(2.0)*SD2) 

* Calculate the angle of diffraction caused by each edge. 

PHIDIFF1=PHI1 

PHID1FF2=PHI2 

* Calculate the change in vertical location of the entering energy bundle, R. 

R1=TAN(PHIDIFF1)*Z 

R2=TAN(PHIDIFF2)*Z 

* Calculate the total distance traveled by the energy bundle before reaching the screen, DIST. 

DIST1=SQRT(Z**2+ Rl**2) 

DIST2=SQRT(Z**2+ R2**2) 

* Calculate the coordinate at which each energy bundle arrives at the observation screen. 

YSCREEN1=Y+R1 

YSCREEN2=Y+R2 

* Increment counter for correct strip on observation screen for ray I . 

IF (ABS(YSCREENl).LE.YSCMAX) THEN 
IF (YSCREEN1.LT.0) THEN 
YDIV =1NT(YSCREEN1/INCR£M)-1 
YDIV=((H/2)+l)+YDIV 
ELSE IF (YSCREEN1.GE.0) THEN 
YDIV =INT(YSCREEN1/INCREM)+1 
YDIV=YDIV+(H/2) 

END IF 

NY(YDIV)=NY(YDIV)+1 
END IF 
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* Increment counter for correct strip on observation screen for ray2. 

IF (ABS(YSCREEN2).LE.YSCMAX) THEN 
IF (YSCREEN2.LT.0) THEN 
YDIV =INT(YSCREEN2/INCR£M)-1 
YDI V=((H/2)+ 1 )+YDIV 
ELSE IF (YSCREEN2.GE.0) THEN 
YDIV =INT(YSCREEN2/INCR£M)+I 
YDlV=YDIV+(H/2) 

END IF 

NY(YDIV)=NY(YDIV)+1 
END IF 

10 CONTINUE 

* After all rays have been fired, store the results in an output file. 

DO 40 1=1, H 

WRITE( 1 5,*)NY(I),(YSCMIN+INCR£M*(I- 1 )),' -\(YSCMIN+INCREM*(I)) 
40 CONTINUE 

END 


* Subroutine that determines x, when erf(x) Is known 
FUNCTION XVAL(ERFX) 

* This subprogram is used to solve for x, knowing erf(x), which is sent from the main program. A 

* truncated infinite series is used to approximate the error function, which can be solved for x 

* using the bisection method. 

* Initialize all variables 
INTEGER NUM,I,NEG 

DOUBLEPRECISION SD,A,B,TOL,P,PO,PI,CHK,ERFX,FP,FA,FB 
REAL XV AL 

DATA TOL/O.OOOl/, PI/3. 141593/, NUM/1000/ 

IF (ERFX.LT.O) THEN 
NEG=1 

ERFX=ABS(ERFX) 

ELSE IF (ERFX.GT.O) THEN 
NEG=2 
END IF 

IF (ERFX.LE.0.992869) THEN 

* Initialize endpoints (x=a, lower endpoint; x=b, upper endpoint) 

A=-l 

B=2.0 

* Initialize counter 
1=0 

* Begin loop to continue for a maximum of NUM iterations 
DO WHILE (I.LT.NUM) 

CHK= ABS((B-A)/2) 
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* Determine midpoint 
P= A + ((B-A)/2) 


* Solve error function for current midpoint and endpoints using a truncated (15 terms) infinite series 

* approximating erf(x). 

FP=(2/SQRT(PI))*(P-(P**3)/(3*l)+(P**5)/(5*2)-(P**7)/(7*6) 
*+(P**9)/(9*24)-(P**ll)/(ll*120)+(P**13)/(13*720)- 
*(P**15)/( 15*5040) +(P**17)/(17*40320)-(P**19)/( 19*362880) 
*+(P**21)/(21*3628800)-(P**23)/(23*39916800))-ERFX 

FA=(2/SQRT(PI))*(A-(A**3)/(3*l)+(A**5)/(5*2)-(A**7)/(7*6) 

*+(A**9)/(9*24)-(A** 1 1 )/( 1 1 * 120)+(A** 13)/(13*720)- 
*(A**15)/(15*5040) + (A**17)/(17*40320) -(A** 1 9)/( 1 9*362880) 
*+(A**21)/(21*3628800)-(A**23)/(23*39916800))-ERFX 

FB=(2/SQRT(PI))*(B-(B**3)/(3* l)+(B**5)/(5*2)-(B**7)/(7*6) 

*+(B**9)/(9*24)-(B**l 1)/(1 1*120)+(B** 13)/( 13*720)- 
*(B**15)/(15*5040) + (B** 1 7)/( 1 7*40320) -(B**19)/(19*362880) 
*+(B**21)/(21*3628800)-(B**23)/(23*39916800))-ERFX 

IF (FP.EQ.0) THEN 
XVAL=P 
GO TO 10 

ELSE IF (CHK.LT.TOL) THEN 
XVAL=(A+B)/2 
GO TO 10 

ELSE IF ((FA*FP).LT.0) THEN 
B=P 

ELSE IF ((FB*FP).LT.0) THEN 
A=P 
END IF 

END DO 

ELSE 
XVAL=2 
END IF 

10 IF (NEG.EQ.l) THEN 
XVAL=-XVAL 
ERFX=-ERFX 
END IF 

END 

* Subroutine that generates Random Numbers 
SUBROUTINE RMARIN(ij, kl) 


C This is the initialization routine for the random number generator ranmar(). NOTE: The seed variables 
C can have values between: 

C 0 <= IJ <= 31328 
C 0<=KL<= 30081 


C The random number sequences created by these two seeds are of sufficient length to complete an entire 
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C calculation with. For example, if sveral different groups are working on different parts of the same 
C calculation, each group could be assigned its own IJ seed. This would leave each group with 30000 
C choices for the second seed. That is to say, this random number generator can create 900 million different 
C subsequences - with each subsequence having a length of approximately 10 A 30. 

C Use IJ = 1802 & KL = 9373 to test the random number generator. The subroutine ranmar should be used 
C to generate 20000 random numbers. Then display the next six random numbers generated multiplied by 
C 4096*4096. If the random number generator is working properly, the random numbers should be* 

C 6533892.0 14220222.0 7275067.0 

C 6172232.0 8354498.0 10633180.0 

C implicit real* 8 (a-h, o-z) 
real*8 u(97), c, cd, cm, s, t 
integer ii, i, j, ij, jj, k, kl, 1, m, i97, j97 
logical test 

common /rasetl/ u, c, cd, cm, i97, j97, test 
test = .false, 
c 

if( IJ .It. 0 .or. IJ gt. 31328 .or. 

S KL .It. 0 .or. KL.gt. 30081 ) then 
write (*, *) ' The first random number seed must have a' 
write (*, *) ' value between 0 and 3 1328/ 
write (*, *) 

write (*, *) ' The second seed must have.a~value between 0* 
write (*, *) ' and 30081/ 
write (*, *) 'Stopping../ 
stop 
endif 
c 

i = mod(IJ/177, 177) + 2 
j = mod(IJ ,177) + 2 
k = mod(KL/169, 178)+ 1 
I = mod(kI, 169) 
c 

do 2 ii = 1,97 
s = 0.0 
t * 0.5 

do 3 jj = 1 , 24 

m - mod(mod(i*j, I79)*k, 179) 

i=j 
j-k 
k = m 

l = mod(53*l+l, 169) 
if (mod(l*m, 64) .ge. 32) then 
s = s + 1 
endif 
t = 0.5 * t 
3 continue 
u(ii) = s 
2 continue 

c= 362436.0/ 16777216.0 
cd = 7654321.0/ 16777216.0 
cm = 16777213.0/16777216.0 

i97 = 97 
J97 = 33 
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c 

test = .true, 
c 

return 

end 

real* 8 function ranmar() 
c 

c This is the random number generator proposed by George Marsaglia 
c in Florida State University Report: FSU-SCRJ-87-50 
c 

c implicit real*8 (a-h, o-z) 
real* 8 u(97), uni, c, cd, cm 
integer i97, j97 
logical test 

common /raset 1/ u, c, cd, cm, i97, j97, test 
c 

if(.not.test) then 
write (*, *) 

write (*, *) 'ranmar error # 1 : must call the’ 
write (*, *) 'initialization routine rmarin before' 
write (*, *) 'calling ranmar.' 
write (*, *) ’Stopping...’ 
stop 
endif 
c 

uni = u(i97) - u(j97) 

if( uni .It. 0.0 ) uni = uni + 1 .0 

u(i97) - uni 

i97 = i97 - l 

if(i97 .eq. 0) i97 = 97 

j97 = j97 - 1 

if(j97 .eq. 0)j97 = 97 

c = c - cd 

if( c .It. 0.0 ) c = c + cm 
uni = uni - c 

if( uni .It. 0.0 ) uni = uni + 1 .0 
c 

ranmar = uni 
c 

return 

end 
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This portion of a program , written in FORTRAN, demonstrates the application of the statistical approach 
to the modeling of diffraction. It predicts the spatial distribution of diffracted energy as it passes through a 
single infinite slit, and arrives at an observation screen some distance away. The subroutines called can 
be found with the first code (for the splitting rays approach) in Appendix A. This particular version is for 

the summing angles approach. 

PROGRAM DIFFRACT2 

C$ NOEXTENSIONS NOWARNINGS 

* Initialize variables used in program 

* Initialize arrays 
DIMENSION NY(1000) 

* Note that array NY is of size H, where H is the number of strips into which observation screen is 

* divided. 

DOUBLE PRECISION LAMDA, W,Z,Y,R1 ,YO,DEL 1 ,DEL2,RANMAR, 

* PI,ERFX 1 ,ERFX2,SD l ,SD2,K t PHI 1 ,PHI2,PHIDIFF, YSCREEN,ARG 1 , 

*YSCMrN,YSCMAX,INCREM,NY,DIST,R,ARG2 

REAL XVAL 

INTEGER N,YDIV,H, NUMRAYS 

DATA LAMDA, W,Z, PI/100.0,60.0,60.25, 3T41593/ 

H=1000 

YO=-W/2 

YSCMAX=80.0 

YSCMIN=-80.0 

INCREM=(YSCMAX-YSCMIN)/H 

NUMRAYS=2000000 

K=(2* PI)/LAMDA 

* Open file that will store output from execution of code. 

OPEN( 1 5,FILE='sumr l OO’^TATUS-'OLD') 

* Initialize the values in the matrix containing the number of energy bundles arriving at screen, which is 

* divided into H strips. 

DO 5 1=1, H 
NY(I)=0 
5 CONTINUE 

* Assign seeds for random number generator. 

1=1 

J= 3 

CALL RMARIN(I,J) 

* Begin Monte-Carlo solution. This loop causes NUMRAYS rays to be fired in randomly and uniformly 

* from aperture. 

DO 10 J=l, NUMRAYS 

* Calculate point of entry of the current energy bundle. 

Y=YO + RANMAR0*W 

* Calculate distance from point of entry of the current energy bundle and the two aperture edges. 
DELl=(W/2)+ABS(Y) 
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DEL2=(W/2)-ABS(Y) 

* Calculate the standard deviation for distribution of diffraction angles. 

SD I =ATAN( 1 0/(2.0*DEL 1 *K)) 

SD2=ATAN( 1 .0/(2. 0* DEL2* K)) 

* Calculate the angle of diffraction caused by each of the aperture edges. The subprogram XVAL must 

* be called in order to determine the value of x, where erf(x) is known. 

ERfXl = 2.0*RANMAR()-1.0 
ARGl=XVAL(ERf XI) 

PHI 1 = ATAN( ARG 1 * SQRT(2 .0)* SD 1 ) 

* Calculate the angle of diffraction caused by the other aperture edge. 

ERFX2= 2.0*RANMAR()-1.0 
ARG2=XVAL(ERFX2) 

PHI2=ATAN(ARG2*SQRT(2.0)*SD2) 

* Calculate the angle of diffraction due to the presence of both sides of the aperture by summing the two 

* angles. 

PHIDIFF=PHI1+PHI2 

* Calculate the change in vertical location of the entering energy bundle, R. 

R=TAN(PHIDIFF)*Z 

* Calculate the total distance traveled by the energy bundle before reaching the screen, DIST. 
DIST=SQRT(Z**2+ R**2) 

* Calculate the coordinate at which each energy bundle arrives at the observation screen. 
YSCREEN=Y+R 

* Increment counter for correct strip on observation screen. 


IF (ABS(YSCREEN).LE.YSCMAX) THEN 
IF (YSCREEN.LT.O) THEN 
YDIV =INT(YSCREEN/TNCREM)-1 
YDlV=((H/2)+l)+YDlV 
ELSE IF (YSCREEN.GE.O) THEN 
YDIV =INT(YSCREEN/INCR£M)+1 
YDIV=YDIV+(H/2) 

END IF 

NY(YDIV)=NY(YDIV)+ 1 
END IF 


10 CONTINUE 

* After all rays have been fired, store the results in an output file. 

DO 40 1=1, H 

WRITE(15,*)NY(I),(YSCMIN+INCREM*(I-1)),' -\(YSCMIN+INCREM*(I)) 
40 CONTINUE 
END 
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This portion of a program, written in FORTRAN, demonstrates the application of the statistical approach 
to the modeling of diffraction . It predicts the spatial distribution of diffracted energy as it passes through a 
single infinite slit, and arrives at an observation screen some distance away . The subroutines called can 
be found with the first code (for splitting rays approach) in Appendix A.. This particular version is for the 
closest edge effect approach . 

PROGRAM DIFFRACT3 

C$ NOEXTENSIONS NOWARNINGS 


* Initialize arrays, note that array NY is of size H, where H is the number of strips into which 

* observation screen is divided. 

DIMENSION NY(IOOO) 

* Initialize variables used in program 

DOUBLE PRECISION LAMDA, W,Z, Y,R 1 , YO,DEL 1 ,DEL2,RANMAR,ARG2 
* PLERFX 1 ,ERFX2,SD 1 ,SD2,K,PHI 1 ,PHI2,PHIDIFF 1 , YSCREEN I ,ARG 1 , 
*YSCMIN,YSCMAX,INCREM,NY,DISTl,PHIDIFF2 t R2,DIST2,YSCREEN2 
REAL XVAL 

INTEGER N,YDIV,H, NUMRAYS 

DATA LAMDA, W,Z, PI/100. 0,60.0, 60. 2S,3rl4 1 593/ 

H=1000 

YO=-W/2 

YSCMAX=80.0 

YSCMIN=~80.0 

INCREM=(YSCMAX-YSCMIN)/H 
NUMRAYS=2 000000 
K=(2* PI)/LAMD A 

* Open file that will store output from execution of code. 

OPEN( 1 5,FILE=’ns 1 00’,STATUS='OLD') 

* Initialize the values in the matrix containing the number of energy bundles arriving at screen, which is 

* divided into H strips. 

DO 5 1=1, H 
NY(I)=0 
5 CONTINUE 

* Assign seeds for random number generator. 

1=1 

J=3 

CALL RMARIN(LJ) 

* Begin Monte-Carlo solution. This loop causes NUMRAYS rays to be fired in randomly and uniformly 

* from aperture. 

DO 10 J= l, NUMRAYS 

* Calculate point of entry of the current energy bundle. 

Y=YO + RANMAR()*W 

* Calculate distance from point of entry of the current energy bundle and the two aperture edges. 
DELl=(W/2)+ABS(Y) 
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DEL2=(W/2)-ABS(Y) 

IF (DEL1.LT.DEL2)THEN 
SD 1 =ATAN( 1 .0/(2. 0* DEL 1 * K)) 

ERFX1= 2.0* RANMAR()- 1 .0 
ARGI=XVAL(ERFX1) 
PHII=ATAN(ARG1 *SQRT(2.0)*SD1) 
PHIDIFF1=PHI1 
R1=TAN(PHIDIFF1)*Z 
DIST1=SQRT(Z**2+ Rl**2) 

YSCREEN 1=Y+R l 

IF (ABS(YSCREENl).LE.YSCMAX) THEN 
IF (YSCREEN 1 .LT.O) THEN 
YDIV =INT(YSCREEN 1/INCREM)- 1 
YDI V=((H/2)+ 1 )+ YDIV 
ELSE IF (YSCREEN l.GE.O) THEN 
YDIV =INT(YSCREENl/INCREM)+l 
YDlV=YDIV+(H/2) 

END IF 

NY(YDIV)=NY(YDIV)+1 
END IF 


ELSE 

SD2=ATAN( 1 .0/(2 .0* DEL2* K)) 

ERFX2= 2.0*RANMAR()-1.0 
ARG2=XVAL(ERFX2) 

PHI2=ATAN(ARG2*SQRT(2.0)*SD2) 

PHIDIFF2=PHI2 
R2=TAN(PHIDIFF2)*Z 
DIST2=SQRT(Z**2+ R2**2) 

YSCREEN2=Y+R2 

IF (ABS(YSCREEN2).LE.YSCMAX) THEN 
IF (YSCREEN2.LT.0) THEN 
YDIV =INT(YSCREEN2/INCREM)-1 
YDIV=((H/2)+ 1 )+YDIV 
ELSE IF (YSCREEN2.GE.0) THEN 
YDIV =INT(YSCREEN2/INCREM)+1 
YDIV=YDIV+(H/2) 

END IF 

NY(YDIV)=NY(YDIV)+ 1 
END IF 
END IF 

10 CONTINUE 

* After all rays have been fired, store the results in an output file. 

DO 40 1=1, H 

WRITE( 1 5,*)NY(I),(YSCMIN+(NCREM*(I- 1 )),' -',(YSCMIN+rNCREM*(I)) 
40 CONTINUE 
END 
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This portion of a program, written in FORTRAN, demonstrates the application of the statistical approach 
to the modeling of diffraction. It predicts the spatial distribution of diffracted energy as it passes through a 
single infinite slit, and arrives at an observation screen some distance away. The subroutines called can 
be found with the first code (for splitting rays approach) in Appendix A. This particular version is for the 
furthest edge effect approach. 

PROGRAM DIFFRACT4 

C$ NOEXTENSIONS NOWARNINGS 


* Initialize arrays, note that array NY is of size H, where H is the number of strips into which 

* observation screen is divided. 

DIMENSION NY(1000) 

* Initialize variables used in program 

DOUBLE PRECISION LAMDA, W,Z, Y,R 1 , YO,DEL 1 ,DEL2,RANMAR, 

* PI.ERFX 1 ,ERFX2,SD 1 ,SD2,K,PHI 1 ,PHI2,PHIDIFF 1 , YSCREEN 1 ,ARG 1 ,ARG2 
*YSCMIN,YSCMAX,INCREM,NY,DIST1,PHIDIFF2,R2,DIST2,YSCREEN2 
REAL XVAL 

INTEGER N,YDIV,H,NUMRAYS 

DATA LAMDA,W,Z,Pl/100.0, 60.0, 60.2*3:141593/ 

H= 1 000 
YO=-W/2 
YSCMAX=80.0 
YSCMIN=-80.0 

INCREM=(YSCMAX-YSCMIN)/H 

NUMRAYS=2000000 

K=(2* PI)/L AMD A 

* Open file that will store output from execution of code. 

OPEN( 1 5,FILE='fs 1 00\STATUS=’OLD’) 

* Initialize the values in the matrix containing the number of energy bundles arriving at screen, which is 

* divided into H strips. 

DO 5 1=1, H 
NY(I)=0 
5 CONTINUE 

* Assign seeds for random number generator. 

1=1 

J=3 

CALL RMARIN(I,J) 

* Begin Monte-Carlo solution. This loop causes NUMRAYS rays to be fired in randomly and uniformly 

* from aperture. 

DO 10 J=l, NUMRAYS 

* Calculate point of entry of the current energy bundle. 

Y=YO + RANMAR0* W 

* Calculate distance from point of entry of the current energy bundle and the two aperture edges. 

DEL 1 =(W/2)+ABS(Y) 
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DEL2=(W/2)-ABS(Y) 

IF (DEL1.GT.DEL2)THEN 
SD 1 =ATAN( 1 ,0/(2.0* DEL 1 * K)) 

ERFX1= 2.0*RANMAR()-1.0 
ARG I =XVAL(ERFX 1 ) 

PHI 1 =ATAN( ARG 1 *SQRT(2.0)*SD 1) 
PHIDIFF1=PHI 1 
Rt=TAN(PHIDIFFl)*Z 
DIST1=SQRT(Z**2+ Rl**2) 
YSCREEN1=Y+R1 

IF (ABS(YSCREENl).LE.YSCMAX) THEN 
IF (YSCREEN 1 .LT.O) THEN 
YDIV =INT(YSCREENl/INCREM)-l 
YDI V=((H/2)+ 1 )+YDIV 
ELSE IF ( YSCREEN 1.GE.0) THEN 
YDIV =INT(YSCREEN1/INCREM)+1 
YDIV=YDIV+(H/2) 

END IF 

NY(YDIV)=NY(YDIV)+ 1 
END IF 


ELSE ' ' 

SD2=ATAN( 1 ,0/(2.0*DEL2*K)) 

ERFX2= 2.0*RANMAR()-1.0 
ARG2=XVAL(ERFX2) 

PHI2=ATAN(ARG2*SQRT(2.0)*SD2) 

PHIDIFF2=PHI2 
R2=TAN(PHIDIFF2)*Z 
DIST2=SQRT(Z* *2+ R2**2) 

YSCREEN2=Y+R2 

IF (ABS(YSCREEN2).LE.YSCMAX) THEN 
IF (YSCREEN2.LT.0) THEN 
YDIV =INT(YSCREEN2/INCR£M)-1 
YDI V=((H/2)+l)+ YDIV 
ELSE IF (YSCREEN2.GE.0) THEN 
YDIV =INT(YSCREEN2/INCREM)+1 
YDIV=YDIV+(H/2) 

END IF 

N Y( YDI V)=N Y (YDI V)+ 1 
END IF 
END IF 

10 CONTINUE 

* After all rays have been fired, store the results in an output file. 

DO 40 1=1 ,H 

WRITE( 1 5,*)NY(I),(YSCMIN+INCREM*(I- 1 )),’ -\(YSCMIN+INCREM*(I)) 
40 CONTINUE 
END 



This FORTRAN program is used to model diffraction through an infinite slit and is based on modified Huy gen- Fresnel 
principle. The phase of an energy bundle upon its arrival to the observation screen is proportional to its optical 
length, and is summed with the phase of the others arriving at the same "bin” of the observation screen. The final sum 
of these phases at a is squared to provide the normalized intensity at that bin. Plotting this intensity of all bins yields a 
diffraction pattern. The subroutine rmarin is not provided here , as it can be found in Appendix A. The following 
defines the variables used in this program. 

* NY is an array that divides the observation screen into strips. Each element in this array serves as a "bin” which 
keeps track of the number of bundles reaching that strip. This array is never actually used in the determination of 
the intensity pattern, but serves to verify that ail of the bins are receiving the same number of energy bundles. 

* LAMDA is the wavelength of the entering monochromatic radiation. 

* W is the width of the rectangular slit. 

* Z is the distance from the aperture to the observation screen. 

* R is the change in vertical location of the entering energy bundle. 

* DIST is the net distance traveled by energy bundle from aperture before being 
intercepted by observation screen. 

* YO is the coordinate of the left hand bottom comer of the rectangular 
aperture. 

* Y is the coordinate of the entering energy bundle in the plane containing the 
aperture. 

* YSCREEN is the coordinate of the energy bundle when it strikes the observation 

screen. Note that the YSCREEN coordinate is such that YSCREEN=0 occurs at the 
aperture center. _ _ 

* YSCMIN, YSCMAX determine the minimum and maximum values of Y for which the 
"number of energy bundles striking” at the observation screen will be 

recorded. 

* H is the number of strips into which the observation screen will be divided. 

* INCREM is the width of each of the strips on the observation screen. 

* RANMAR calls a random number from the subroutine rmarin(ij, kl). 

* K is the wavenumber, given by 2*PI/LAMDA. 

* YDIV is the variable used to determine the "bin” into which to store 
energy bundles incident to the observation screen. 


PROGRAM SLITMONTECARLO 
dimension ny( 10000), beta( 10000) 

double precision lamda, w, z, y, r, yo, ranmar, pi, phidiff, yscreen, beta, k, yscmax, increm, dist, 
41 phimin, phimax, delphi,betaold,yscmin,obf 
integer n,ydiv,h,numrays,q,I j 
open(0 1, file-out.dat', status- old’) 

data lamda, w,z,pi/0.00000058, 0.0003,16.0,2.0,3.14159/ 
yo=-w/2.0d0 

* Note that h must be evenly divisible by 2 
h= 10000 

* Limit range of interest; specifying min/max y of interest 
yscmax=0.009d0 

yscmin=-yscmax 

numrays=l00000 

k=2.0d0*pi/lamda 
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* Initialize the values in the matrix containing the number of energy bundles arriving at screen in a certain strip, 
do 3 j=l,h 

ny(j)=0.0d0 
beta(j)=0.0d0 
3 continue 

* Initialize random number generator. 
i=l 

j=3 

increm=(yscmax-yscmin)/h 
call rmarin(ij) 


do 5 q=l,numrays 

* Calculate the random point of entry of current energy bundle. Note that the following restricts the direction of 

* emission so that ail energy bundles will arrive at the observation screen, but maintains a diffuse emission where all 

* diffraction angles within the limited range will be equally probable. 

y=yo+ranmar()*w 

phimax=datan((yscmax)/z) 

phidiff=-phimax+2*phimax*ranmar() 

r=dtan(phidiff)*z 

yscreen=y+r 

dist=dsqrt(z* * 2+r* *2) 

* Increment the number falling into a bin, keeping with the numbering scheme. 
if(abs(yscreen).le.yscmax)then 

if(yscreen.lt,O.OdO)then 
ydiv=int(abs(yscreen/increm))+ 1 
ydiv=(h/2)-ydiv+l 
else if (yscreen.ge.O.OdO) then 
ydiv=int(yscreen/increm)+ 1 
ydiv=ydiv+(h/2) 
end if 

ny(ydiv)=ny(ydiv)+ 1 

* Increment the phase of the appropriate bin. 
betaold=beta(ydiv) 

* Define the obliquity factor. 
obf=0.5d0*(l .OdO+dcos(phidiff)) 

* Note that here dcos(2d0*dist*pi/lamda)) is used instead of dsin (2d0*dist*pi/lamda)), as defined by equation 4.33. 

* Either provide similar results, whereby the peaks fall within the envelope of the analytical solution, however use of 

* cos results in a pattern with the central oscillation at a peak, where use of sin results in the central oscillation at a 

* valley. 

beta(ydiv)=(obf*dcos(2d0*dist*pi/lamda))/(2d0*dist*pi/lamda) 

beta(ydiv)=beta(ydiv)+betaold 

end if 

5 continue 
do 4 j=l,h 

write(01,*)sngl(atan((yscmin+increm^j)/z)),sngl(abs(beta0)**2)) 

4 continue 
end 



/*The following code was written in C programming language. It serves to define the geometry for the 
radiative model of the current CERES telescope and can be modified to model potential next-generation 
instruments. It calls C +> library functions written by TRG doctoral student, Felix Nevarez. This code 
generates output which can be used to determine the Optical Point Spread Function of the CERES 
instrument. Each surface within the telescope is modeled, and is referenced by a combination of letters and 
numbers. This labeling corresponds to the labeled Figure 5. 1 (a) of the CERES telescope found in this 
thesis. Note that all dimensions are in inches, scaled by a factor of ten ( 1 0). */ 

^include "raytracer.h" 

^include "math.h" 

^define PI 3.14159265359 

main() 

{ 

/*s is the scaling factor to scale all telescope dimensions 41 / 
double s= 10; 

/*sh is the shift added by the insertion of shims to achieve best focus 41 / 
double sh=. 01 10 ; 

/*apw is the width of a precision aperture 41 / 
double apw=0.0296; 

/* d is the distance between precision apertures*/ 
double d=0,01; 

/*theta is direction 1 of incoming rays*/ * ' 
double theta; 

/*phi is direction 2 of incoming rays*/ 
double phi; 

/* Begin defining the CERES telescope geometry*/ 

/*CY1 (Baffle)*/ 

makeACylinder(0.56l*s); 

setProperties(0.9,0.1); 

setClipPlane (0.0,0.0,(0.0+o),0.0,0.0, 1 .0); 

setClipPlane(0.0,0.0,( 1 .77 1 +o)*s,0.0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 

/*R1 (Baffle)*/ 

makeARing(0.0, 0. 0, (0.0+o)*s, 0.0, 0.0, 1.0,0. 4665*s, 0. 56 10*s); 

setProperties(0.99,0. 1); 

noOutput(); 

nextOBJ(); 

/*R2(Baffle)*/ 

makeARing(0.0, 0.0, (0.39+o)*s ,0.0,0.0,1.0,0.4570*s, 0.56 10*s); 

setProperties(0.99,0. 1 ); 

noOutput(); 

nextOBJ(); 

/*R3(Baffle)*/ 

makeARing(0.0,0.0,(0.770+o)*s,0.0,0.0, 1 .0,0.4475*5,0.56 1 0*s); 

setProperties(0.99,0. 1 ); 

noOutput(); 

nextOBJ(); 


/*R4(Baffle)*/ 
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makeARing(0.0,0.0,( l . 149+o)*s ,0.0, 0.0,1 . 0,0.4385^,0.56 10*s); 

setProperties(0.99,0. 1); 

noOutputO; 

nextOBJ(); 

/*R5(Baffle)*/ 

makeARing(0.0,0.0,( 1 .529+o)*s ,0.0, 0.0, l .0,0.4290*s,0.56 1 0*s); 

setProperties(0.99,0. 1 ); 

noOutput(); 

nextOBJ(); 

/* Cl (Baffle)*/ 

makeACone(0.3239*s,0.56 1 *s); 

setProperties(0.99,0. 10); 

set0rigin(0.0,0.0,(2.0949+o)*s); 

setClipPlane(0. 0, 0. 0,(1. 77 l+o)*s, 0.0,0. 0,1.0); 

setClipPlane(0.0,0.0,( 1 .8220+o)*s,0.0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 

/*CY2(BaffIe)*/ 

makeACylinder(0.472*s); 

setProperties(0.99,0. 1); - ' 

setClipP!ane(0.0,0.0,( 1 .8220+o)*s,0.0,0.0, 1 .0); 

setClipPlane(0.0,0.0,( 1 .8540+o)*s,0.0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 

/*R6(Baffle)*/ 

makeARing(0.0,0.0,(1.854+o)*s,0.0,0.0,1.0,0.3995*s,0.472*s); 

setProperties(0.99,0. 1); 

noOutputO; 

nextOBJ(); 

/*C2(CapReflector)*/ 
makeACone(0.952*s,0.424*s); 
setProperties(0.99,0. 1); 
set0rigin(0.0,0.0,(0,9567+o)*s); 
setClipPlane(0. 0,0. 0,(1 .854+o)*s,0. 0,0.0, 1.0); 
setClipPlane(0.0,0.0,( 1 .909+o)*s,0. 0,0.0,- 1 .0); 
noOutputO; 
nextOBJ(); 

/*C3(CapReflector)*/ 
makeACone(0. 1904*s,0.424*s); 
setProperties(0.00 1 ,0.9); 
set0rigin(0.0,0.0,(2.0994+o)*s); 
setClipPlane(0.0, 0.0, (l.909+o)*s, 0.0,0. 0,1.0); 
setClipPlane(0.0,0.0,(l.920+o)*s,0. 0,0.0, -1.0); 
noOutputO; 
nextOBJ(); 

/*C4(CapReflector)*/ 
makeACone(0.952*s,0.424*s); 
setProperties(0.99,0. 1 ); 
set0rigin(0.0,0.0,( 1 .0227+o)*s); 
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setClipPlane(0.0,0.0,( 1 .92+o)*s,0. 0,0.0, 1 .0); 
setClipPlane(0.0,0.0,( 1 .975+o)*s,0.0,0.0,- 1 .0); 
noOutput(); 
nextOBJ(); 

/*C5(CapReflector)*/ 
makeACone(0. I904*s,0.424*s); 
setProperties(0.00 1,0.9); 
setOrigin(0. 0,0. 0,(2. 1654+o)*s); 
setClipPlane(0.0,0.0,( 1 .975+o)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0,0.0,( 1 .986+o)*s,0.0,0.0,- 1 .0); 
noOutput(); 
nextOBJ(); 

/*C6(CapReflector)*/ 
makeACone(0.952*s,0.424*s); 
setProperties(0.99,0. 1 ); 
set0rigin(0.0,0.0,( 1 .0887+o)*s); 
setClipPlane(0.0,0.0,( 1 .986+o)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0,0.0,(2.041+o)*s,0.0,0.0,- 1 .0); 
noOutputQ; 
nextOBJ(); 

/*C7(CapReflector)*/ 
makeACone(0. 1904*s, 0.424*$); 
setProperties(0.00 1 ,0.9); 
setOrigin(0. 0,0. 0,(2. 23 14+0)*$); 
setClipP!ane(0. 0,0.0, (2.04 1 +o)*s,0.0,0.0,- 1 .0); 
setClipPlane(0.0,0.0,(2.052+o)*s,0.0,0.0, 1 .0); 
noOutputO; 
nextOBJ(); 

/*C8(CapReflector)*/ 
makeACone(0.952*s,0.424*s); 
setProperties(0.99,0. 1); 
set0rigin(0.0,0.0,( 1 . 1 547+o)*s); 
setClipPlane(0. 0,0. 0,(2. 052+o)*s, 0.0,0. 0,1.0); 
setClipPlane(0. 0,0. 0,(2. 1 07+o)*s,0. 0,0.0,- 1 .0); 
noOutput(); 
nextOBJ(); 

/*C9(CapReflector)*/ 
makeACone(0. 1904*s,0.424*s); 
setProperties(0.00 1 ,0.9); 
set0rigin(0.0,0.0,(2.2974+o)*s); 
setClipPIane(0. 0,0. 0,(2. 1 07+o)*s,0.0,0.0, 1 .0); 
setClipPIane(0. 0,0.0, (2.1 l8+o)*s,0.0,0.0,-1.0); 
noOutput(); 
nextOBJ(); 

/* C l 0(Cap Reflector)* / 
makeACone(0.952*s,0.424*s); 
setProperties(0.99,0. 1 ); 
set0rigin(0.0,0.0,(1.2207+o)*s); 
setClipP!ane(0.0,0.0,(2. 1 1 8+o)*s,0.0,0.0, l .0); 
setClipP!ane(0.0,0.0,(2. 173+o)*s,0.0,0.0,-1.0); 
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noOutput(); 

nextOBJ(); 

/♦Cll(CapReflector)*/ 
makeACone(0. 1904*s,0.424*s); 
setProperties(0.001 ,0.9); 
set0rigin(0.0,0.0,(2.363+o)*s); 
setClipPlane(0. 0,0. 0,(2. l73+o)*s,0.0,0.0, l .0); 
setClipPlane(0.0,0.0,(2. 1 84+o)*s,0.0,0.0,- 1 .0); 
noOutput(); 
nextOBJ(); 

/*C 1 2(CapReflector)*/ 
makeACone(0.924*s,0.41 l*s); 
setProperties(0.99,0. 1 ); 
set0rigin(0.0,0.0,( 1 .287+o)*s); 
setClipPlane(0.0,0.0,(2. 1 84+o)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0,0.0,(2.2 1 l+o)*s,0.0,0.0,- 1 .0); 
noOutput(); 
nextOBJQ; 

/*R7(CapReflector)*/ 

makeARing(0.0,0.0,(2.2 1 1 +o)*s,0.0,0.0, l .0,0,4 1 1* s,0.4745*s); 

setProperties(0.99,0. 1); 

noOutputO; 

nextOBJ(); 

/*CY3(Cap Reflector)*/ 

makeACylinder(0.4745*s); 

setProperties(0.99,0. 1 ); 

setClipPlane(0.0,0.0,(2.2 1 1 +o)*s,0.0,0.0, 1 .0); 

setClipPIane(0.0,0.0,(2.221+o)*s,0.0,0.0,-l.0); 

noOutputO; 

nextOBJ(); 

/*R8(CapReflector)*/ 

makeARing(0.0,0.0,(2221+o)*s,0.0,0.0 f 1.0,0.4745*s,0.522*s); 

setProperties(0.0 1 ,0.9); 

noOutput(); 

nextOBJ(); 

/*CY4(CapReflector)*/ 

makeACylinder(0.522*s); 

setProperties(0. 01,0.9); 

setClipPIane (0.0, 0.0, (2.22 1 +o)*s,0.0,0.0, 1 .0); 

setClipPlane(0.0,0.0,(2.263+o)*s,0. 0,0.0,- 1 .0); 

noOutputO; 

nextOBJ(); 

/*R9(CapReflector)*/ 

makeARing(0.0,0.0,(2.263+o)*s , 0.0, 0. 0,1. 0,0. 4725 *s,0.5220*s); 

setProperties(0.0 1 ,0.9); 

noOutputO; 

nextOBJQ; 


/*CY5(TelescopeHousing)*/ 
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makeACylinder(0.4725*s); 

setProperties(0.99,0. 1 ); 

setClipPlane (0.0,0.0,(2.263+o)*s,0.0,0.0,1.0); 

setClipPlane(0. 0,0. 0,(2. 386+o)*s, 0.0, 0.0,- 1 .0); 

noOutput(); 

nextOBJQ; 

/*D1 (Spider)*/ 

makeADisk(0.0, 0.0,(2.263+o)*s, 0.0,0. 0,1.0,0.185*s); 

setProperties(0.99,0. 1 ); 

noOutputQ; 

nextOBJ(); 

/♦legl(Spider)*/ 

makeAPlane(0. 0,0. 0,(2. 263+0. 1 13+o)*s,0.0,0.0, 1 .0); 
setProperties(0.99,0. 1 ); 

setClipPlane(0.03*s,-0.4725*s,(2.263+0.l 13+o)*s,0.0,1.0,0.0); 
setClipPlane(0.03*s, -0.4725*s, (2.263+0. 1 13+o)*s,-l.0, 0.0, 0.0); 
setClipPlane(-0.03*s,0. 0,(2. 263+0.1 13+o)*s,l .0,0. 0,0.0); 
setClipPlane(0*s,-0. 185*s,(2.263+0. 1 13+o)*s,0. 0,-1 .0,0.0); 
noOutput(); 
nextOBJ(); 

/*leg2(Spider)*/ 

makeAPlane(0.0, 0.0, (2.263+0. 1 1 3+o)*s,0.0,0.0, 1 .0); 
setProperties(0.99,0. 1); 

setClipPlane(0. 1 60*s, 0.0925*s, (2.263+0. 113+o)*s, 0.866, 0.5, 0.0); 
setClipPIane(0.459*s,0.23045*s,(2.263+0.1 13+o)*s, -0.5, 0.866, 0.0); 
setClipPlane(0.459*s,0.23045*s, (2.263+0. 1 1 3+o)*s, -0.866, -0.5, 0.0); 
setClipPlane(0.429*s,0.28225*s, (2.263+0.1 13+o)*s, 0.5, -0.866, 0.0); 
noOutput(); 
nextOBJQ; 

/*leg3(Spider)*/ 

makeAPlane(0.0, 0.0, (2.263+0. 1 1 3+o)*s,0.0,0.0, 1 .0); 
setProperties(0.99,0. 1); 

setCiipPlane(-0.160*s,0.0925*s, (2.263+0. 113+o)*s,-0. 866, 0.5, 0.0); 
setClipPlane(-0.459 <t s,0.23045*s, (2.263+0. 1 1 3+o)*s, 0.5, 0.866, 0.0); 
setClipPlane(-0.459*s,0.23045*s,(2.263+0. 1 13+o)*$, 0.866,-0.5,0.0); 
setClipPlane(-0.429*s,0.28225*s,(2.263+o+0.1 13)*s, -0.5, -0.866, 0.0); 
noOutputO; 
nextOBJ(); 

/* s ide 1 to leg 1 (Spider)* / 

makeAPlane(0.03*s, 0.0, (2.263+0. 113+o)*s, -0.9659, 0.0, 0.2588); 
setProperties(0.00 1 ,0.9); 

setClipPlane(0.03*s, -0.4725*s, (2.263+0. 113+o)*s, 0.0, 1.0, 0.0); 
setClipPlane(0.0, 0.0, (2.263+0. 1 1 3+o)*s,0.0,0.0,- 1 .0); 
setClipPlane(0.0,0.0,(2.263+o)*s,0. 0,0. 0,1.0); 
setClipPlane(0*s,-0. 1 85*s, (2.263+0. 1 l3+o)*s,0.0,- 1 .0,0.0); 
noOutput(); 
nextOBJ(); 

/*side2toleg 1 (Spider)*/ 

makeAPlane(-0.03*s, 0.0, (2.263+0. 113+o)*s, -0.9659, 0.0, -0.2588); 
setProperties(0 .00 1,0.9); 
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setClipPlane(0.03*s, -0.4725*s, (2.263+0. 113+o)*s, 0.0, 1.0, 0.0); 
setClipPlane(0.03*s, 0.0,(2.263+0. 1 13+o)*s,0.0.0.0,- 1 .0); 
setCIipPlane(0. 0,0 .0,(2. 263+o)*s,0. 0,0.0, 1 .0); 
setCIipPlane(0*s,-0. 185*s, (2.263+0. 1 13+o)*s,0. 0,-1. 0,0.0); 
noOutput(); 
nextOBJ(); 


/*sideltoleg2(Spider)*/ 

makeAPtane(0.429*s,0.28225*s, (2.263+0. 113+o)*s, 0.48295,-0.866,0.2588); 
setProperties(0.00 1 ,0.9); 

setClipPlane(0.160*s,0.0925*s, (2.263+0. 113+o)*s, 0.866, 0.5, 0.0); 
setClipPlane(0.0,0.0,(2.263+o)*s,0.0,0.0,1.0); 

setClipPlane(0.4092*s,0.23625*s,(2.263+o+0. 1 1 3)*s,-0. 866,-0.5,0.0); 

setCIipPlane(0. 0,0.0, (2.263+O+0.1 13)*s,0.0,0.0,-l.0); 

noOutput(); 

aextOBJ(); 

/*side2toleg2(Spider)*/ 

makeAPIane(0.459*s,0.23045*s, (2.263+0. 113+o)*s, -0.48295, 0.866,0.2588); 
setProperties(0.00 1,0.9); 

setClipPlane(0.160*s,0.0925*s, (2.263+0. 113+o)*s, 0.866, 0.5, 0.0); 
setClipPlane(0.0,0.0,(2.263+o)*s,0.0,0.0, 1 .0)^ - 

setClipPlane(0.4092*s, 0.23625*s, (2.263+0. 113+o)*s, -0.866, -0.5, 0.0); 

setClipPlane(0.0, 0.0,(2.263+0. 1 1 3+o)*s,0.0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 


/*sideltoleg3(Spider)*/ 

makeAPlane(-0.459*s,0.23045*s, (2.263+0. 1 13+o)*s,0.48295, 0.866, 0.2588); 
setProperties(0.00 1 ,0.9); 

setClipPlane(-0.160*s,0.0925*s, (2.263+0. 113+o)*s, -0.866,0.5,0.0); 

setClipPlane(-0.4092*s,0.23625*s,(2.263+o+0.1 13)*s,0.866,-0.5,0.0); 

setClipPlane(0.0,0.0,(2.263+o)*s,0.0,0.0, 1 .0); 

setClipPlane(0.0,0.0,(2.263+o+0.1 13)*s,0.0,0.0,-1.0); 

noOutput(); 

nextOBJ(); 

/*side2toleg3(Spider)*/ 

makeAPlane(-0.429*s, 0.28225*s, (2.263+0. 113+o)*s, -0.48295, -0.866, 0.2588); 
setProperties(0.00 1 ,0.9); 

setClipPlane(-0.160*s, 0.0925*s, (2.263+0.1 13+o)*s,-0.866, 0.5, 0.0); 

setClipPlane(-0.4092*s,0.23625*s, (2.263+0. 1 1 3+o)*s, 0.866,-0.5,0.0); 

setClipPlane(0.0,0.0,(2.263+o)*s,0.0,0.0, 1 .0); 

setClipPlane(0.0, 0.0,(2.263+0+0. 1 1 3)*s,0.0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 

/*CY6(Spider)*/ 

makeACylinder(0. 1 85* s); 

setProperties(0.99,0. 1); 

setClipPtane (0.0,0.0,(2.263+o)*s,0.0,0.0, 1 .0); 

setClipPlane(0.0,0.0,(2.506+o)*s,0.0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 
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/*CY7(Spider)*/ 

makeACy linder(0. 1 6*s); 

setProperties(0.99,0. 1 ); 

setClipPlane (0.0,0.0,(2.448+o)*s,0.0,0.0,1.0); 

setClipPlane(0.0,0.0,(2.4762+o)*s,0.0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 

/*C13(Spider)*/ 
makeACone(0.2205*s,0. 185*s); 
setProperties(0.99,0. 1 ); 
set0rigin(0.0,0.0,(2.2855+o)*s); 
setClipPlane(0.0,0.0,(2.4762+o)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0, 0.0, (2.506+o)*s,0. 0,0.0, -1.0); 
noOutputO; 
nextOBJ(); 

/*M2(Secondary Mirror)*/ 

makeASphere(1.3 18*s); 

setProperties(0.00 1,1.0); 

setOrigii^O.OAOXl . 14+o)*s); 

setClipPlane(0.0,0.0,(2.448+o)*s,0.0,0.0, 1 .0); - 

noOutputQ; 

nextOBJ(); 

/*C14(TelescopeHousing)*/ 
makeACone(0.2483*s,0.4725*s); 
setProperties(0.99,0. 1 ); 
set0rigin(0.0,0.0,(2.6343+o)*s); 
setClipPIane(0. 0,0.0, (2. 386+o)*s,0. 0,0. 0,1.0); 
setClipP!ane(0.0,0.0,(2.448+o)*s,0.0,0.0,- 1 .0); 
noOutput(); 
nextOBJ(); 

/* R l O(TelescopeHousing)*/ 

makeARing(0.0,0.0,(2.448+o)*s ,0.0, 0.0, 1 .0,0.3545*s,0.4295*s); 

setProperties(0.99,0. 1); 

noOutput(); 

nextOBJ(); 

/*CY8(TelescopeHousing)*/ 

makeACy linder(0.4295*s); 

setProperties(0.99,0. 1 ); 

setClipPlane (0.0,0.0,(2.448+o)*s,0.0,0.0, 1 .0); 

setCIipPlane(0. 0,0.0, (2. 498+o)*s,0. 0,0,0,- 1.0); 

noOutput(); 

nextOBJ(); 

/*R1 l(TelescopeHousing)*/ 

makeARing(0.0,0.0,(2.498+o)*s ,0.0,0.0,1.0,0.365*s,0.4295*s); 

setProperties(0.99,0.1); 

noOutput(); 

nextOBJ(); 


/*C15(TelescopeHousing)*/ 
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makeACone(0.2807*s,0.486*s); 

setProperties(0.99,0. 1 ); 

set0rigin(0.0,0.0,(2.2873+o)*s); 

setClipPlane(0.0,0.0,(2.498+o)*s,0.0,0.0, 1 .0); 

setClipPlane(0. 0,0.0, (2. 568+o)*s,0. 0,0.0,- 1 .0); 

noOutput(); 

nextOBJ(); 

/*C16(TeIescopeHousing)*/ 

makeACone(0.2807*s,0.486*s); 

setProperties(0.99,0. 1); 

set0rigin(0.0,0.0,(2.8487+o)*s); 

setClipPlane(0. 0,0.0, (2. 568+o)*s,0. 0,0.0, 1 .0); 

setClipPlane(0.0,0.0,(2.638+o)*s,0.0,0.0,-1.0); 

noOutput(); 

nextOBJ(); 

i*C 1 7(TelescopeHousing)*/ 
makeACone(0.2807*s,0.486*s); 
setProperties(0.99,0. 1); 
set0rigin(0.0,0.0,(2.4273+o)*s); 
setClipPlane(0.0,0.0,(2.638+o)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0,0.0,(2.708+o)*s,0.0,0.0,- 1 .01; , 
noOutput(); 
nextOBJ(); 

/*C18(TelescopeHousing)*/ 
makeACone(0.2807*s,0.486*s); 
setProperties(0.99,0. 1); 
set0rigin(0.0,0.0,(2.9887+o)*s); 
setClipPlane(0.0,0.0,(2.708+o)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0,0.0,(2.778+o)*s,0.0,0.0,- 1 .0); 
noOutputQ; 
nextOBJ(); 

/*CY9(TelescopeHousing)*/ 

makeACylinder(0.365*s); 

setProperties(0. 99,0.1); 

setClipPlane (0.0,0.0,(2.778+o)*s,0.0,0.0,1.0); 

setClipPIane(0.0,0.0,(2.83 l+o+sh)*s,0.0,0.0,-1.0); 

noOutputQ; 

nextOBJQ; 

/'“Ml (PrimaryMirror)*/ 

makeASphere( 1 .446*s); 

setProperties(0.00 1,1.0); 

set0rigin(0.0,0.0,(l.432+o+sh)*s); 

setClipPlane(0.0, 0.0, (2.83 1 +o+sh)*s,0.0,0.0, 1 .0); 

setClipPlane(0.0,0.0,(2.867+o+sh)*s,0.0,0.0,-l.0); 

noOutputQ; 

nextOBJ(); 

/*C 1 9(PrimaryMirrorInsert)*/ 
makeACone(0.0466*s,0. 174*s); 
setProperties(0.00 1 ,0.9); 
set0rigin(0.0,0.0,(2.820+o+sh)*s); 
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setClipPlane(0.0,0.0,(2.857+o+sh)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0,0.0,(2.867+o+sh)*s,0.0,0.0,- 1 .0); 
noOutput(); 

/*C20(PrimaryMirrorInsert)*/ 
makeACone(0.082*s,0. 135*s); 
setProperties(0.00 1 ,0.9); 
set0rigin(0.0,0.0,(2.939+o+sh)*s); 
setClipPlane(0.0,0.0,(2.857+o+sh)*s,0.0,0.0, 1 .0); 
setCIipPlane(0.0,0.0,(2.915+o+sh)*s,0.0,0.0 f -1.0); 
noOutput(); 

/*C21(PrimaiyMirrorInsert)*/ 
makeACone(0. 109*s,0.0509*s); 
setProperties(0.00 1 ,0.9); 
set0rigin(0.0,0.0,(2.83 l+o+sh)*s); 
setClipPlane(0.0,0.0,(2.9 1 5+o+sh)*s,0.0,0.0, 1 .0); 
setClipPlane(0.0,0.0,(2.9404+o+sh)*s,0.0,0.0,- 1 .0); 
noOutput(); 

/♦Precision Aperture*/ 

/*If number of precision apertures =1*/ 

/*D2 lprecisionaperture*/ 

makeADisk(0.0,0.0,(2.9404+o+sh)*s,0.0,0.0,-l.G,0.3*s); 
setProperties(0.00 1 ,0.9); 
setClipPlane(0.0148*s,0. 0,0. 0,-1. 0,0. 0,0.0); 
setClipPIane(0.0,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPlane(0.0,0.0296*s,0.0, l .0,- 1 .0,0.0); 
setClipPlane(-0.0148*s, 0.0, 0.0,1. 0,0.0, 0.0); 
setClipPlane(0. 0,-0. 0296*s, 0.0, 1 .0, 1 .0,0.0); 
setClipPlane(0.0,-0.0296*s,0.0,-1.0,l .0,0.0); 
clipAsAHole(); 


/*If number of precision apertures =2*/ 

/*D2 2precisionapertures*/ 

makeADisk(0.0,0.0,(2.9404+o+sh)*s,0.0,0.0,- 1 .0,0.03950*s); 
setProperties(0.001,0.9); 

setClipPlane((0.0l48+(apw+d)/2)*s,0. 0,0. 0,-1. 0,0. 0,0.0); 
setClipPlane((apw+d)*s/2,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPIane((apw+d)*s/2,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPlane((((apw+d)/2)-0.0 148)*s,0. 0,0.0, 1 .0,0. 0,0.0); 
setClipPlane((apw+d)*s/2,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane((apw+d)*s/2,-0.0296*s,0.0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

setClipPlane((0.0148-(apw+d)/2)*s,0.0,0. 0,-1. 0,0. 0,0.0); 
setClipPlane(-(apw+d)*s/2,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPIane(-(apw+d)*s/2,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPIane((-((apw+d)/2)-0.0 1 48)*s,0.0,0.0, 1 .0,0.0, 0.0); 
setClipPlane(-(apw+d)*s/2,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane(-(apw+d)*s/2,-0.0296*s,0.0,- 1 .0, 1 .0,0.0); 
clipAsAHoleQ; 

/*If number of precision apertures =3*/ 

/*D2 3precisionapertures*/ 

makeADisk(0. 0,0. 0,(2. 9404+o+sh)*s, 0.0, 0.0, -1.0, 0.3 *s); 
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setProperties(0.00 1 ,0.9); 
setClipPlane(0.0148*s,0.0,0.0,-1.0,0.0,0.0); 
setClipPlane(0.0,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPlane(0.0,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setCIipPlane(-0.0 148*s,0.0,0.0, 1 .0,0. 0,0.0); 
setClipPlane(0. 0,-0. 0296*s, 0.0, 1 .0, 1 .0,0.0); 
setClipP!ane(0. 0,-0. 0296* s, 0.0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

setClipP!ane((0.0148+apw+d)*s,0.0,0.0,-1.0,0.0,0.0); 
setClipPlane((apw+d)*s,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setCiipPlane((apw+d)*s,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setCIipPlane((apw+d-0.0148)*s,0.0,0.0,1.0,0.0,0.0); 
$etClipPlane((apw+d)*s,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane((apw+d)*s,-0.0296*s,0.0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

setClipPiane((0.0 1 48-apw-d)*s,0. 0,0.0,- 1 .0,0. 0,0.0); 
setClipPlane(-(apw+d)*s,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPlane(-(apw+d)*s,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPlane((-apw-d-0.0148)*s,0.0,0.0,1.0,0.0,0.0); 
setClipPlane(-(apw+d)*s,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane(-(apw+d)*s,-0.0296*s,0.0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

/* If number of precision apertures =4*/ 

/*D2 4precisionapertures*/ 

makeADisk(0.0,0.0,(2.9404+o+sh)*s,0.0,0.0,- 1 .0,0. 1 *s); 
setProperties(0.001,0.9); 

setCtipPIane((0.0l48+(apw+d)/2)*s,0. 0,0. 0,-1. 0,0. 0,0.0); 
setClipPlane((apw+d)*s/2,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPlane((ap\v+d)*$/2,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPlane((((apw+d)/2)-0.0148)*s,0.0,0.0, 1 .0,0.0, 0.0); 
setClipPlane((apw+d)*s/2,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane((apw+d)*$/2,-0.0296*s,0.0,-1.0,1.0,0.0); 
clipAsAHole(); 

setClipPlane((0.0148-(apw+d)/2)*s,0. 0,0.0,- 1 .0,0. 0,0.0); 
setClipP!ane(-(apw+d)*s/2,0.0296*s,0.0,-l. 0,-1. 0,0.0); 
setClipP!ane(-(apw+d)*s/2,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPiane((-((apw+d)/2)-0.0148)*s,0, 0,0. 0,1. 0,0. 0,0.0); 
setClipPlane(-(apw+d)*s/2,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane(-(apw+d)*s/2,-0.0296*s,0.0,- 1 .0, l .0,0.0); 
clipAsAHole(); 

setClipPIane((0.0148+3*(apw+d)/2)*s,0.0,0.0,-1.0,0.0,0.0); 
setClipPlane((apw+d)*3*s/2,0.0296*s,0. 0,-1. 0,-1. 0,0.0); 
setClipPlane((apw+d)*3 *s/2,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPlane(((3*(apw+d)/2)-0.0148)*s,0. 0,0. 0,1. 0,0. 0,0.0); 
setClipP!ane((apw+d)*s*3/2,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane((apw+d)*3*s/2,-0.0296*s, 0.0, -t. 0,1, 0,0.0); 
clipA$AHole(); 

setClipPlane((0.0148-3*(apw+d)/2)*s,0. 0,0. 0,-1. 0,0. 0,0.0); 
setClipPlane(-(apw+d)*3*s/2,0.0296*s,0.0,-1.0,-1.0,0.0); 
setCIipPlane(-(apw+d)*3*s/2,0.0296*s,0. 0,1. 0,-1. 0,0.0); 
setCIipPlane((-((apw+d)*3/2)-0.0 148)*s,0.0,0.0, 1 .0,0. 0,0.0); 
setClipPlane(-(apw+d)*s*3/2,-0.0296*s,0. 0,1. 0,1. 0,0.0); 
setClipPlane(-(apw+d)*s*3/2,-0.0296*s,0. 0,-1. 0,1. 0,0.0); 
clipAsAHole(); 
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/*If number of precision apertures =5*/ 

/* D2 Sprecisionapertures*/ 

makeADisk(0.0,0.0,(2.9404+o+sh)*s,0.0,0.0,- 1 .0,0.03950*s); 
setProperties(0.00 1 ,0.9); 
setCIipPIane(0.0148*s, 0.0, 0.0,- 1.0,0. 0,0.0); 
setClipP!ane(0.0,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPIane(0.0,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPlane(-0.0148*s,0.0,0.0,1.0,0.0,0.0); 
setClipPlane(0.0,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane(0. 0,-0. 0296*s, 0.0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

setCtipPlane((0.0148+ap\v-Hi)*s,0.0,0.0,-1.0,0.0,0.0); 
setClipPlane((apw+d)*s,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPlane((apw+d)*s,0.0296*s,0.0, 1 .0,- i .0,0.0); 
setClipPlane((apw+d-0.0 1 48)*s,0. 0,0.0, 1 .0,0. 0,0.0); 
setClipPlane((apw+d)*s,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane((apw+d)*s,-0.0296*s,0,0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

setClipPlane((0.0148-apw-d)*s,0. 0,0. 0,-1. 0,0. 0,0.0); 
setClipPlane(-(apw+d)*s,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPlane(-(apw+d)*s,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPlane((-apw-d-0.0148)*s,0. 0,0. 0,1. 0,0. 0,0.0); 
setClipPlane(-(apw+d)*s,-0.0296*s,0.0, 1 .0, 1 .0,0*.0); 
setClipPlane(-(apw+d)*s,-0.0296*s,0.0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

setClipPIane((0.0148+2*(apw+d))*s,0.0,0.0,-1.0,0.0,0.0); 
setClipPlane((apw+d)*2*s,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setCIipPlane((apw+d)*2*s, 0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPlane((2*(apw+d)-0.0148)*s, 0.0, 0.0,1. 0,0.0, 0.0); 
setCIipPlane((apw+d)*2*s,-0.0296*s,0.0, 1 .0, 1 .0,0.0); 
setClipPlane((apw+d)*2*$,-0.0296*s,0.0,- 1 .0, 1 .0,0.0); 
ciipAsAHole(); 

setClipPlane((0.0148-2*(apw+d))*s,0. 0,0. 0,-1. 0,0. 0,0.0); 
setClipPlane(-(apw+d)*2*s,0.0296*s,0.0,- 1 .0,- 1 .0,0.0); 
setClipPlane(-(apw+d)*2*s,0.0296*s,0.0, 1 .0,- 1 .0,0.0); 
setClipPIane((2*(-apw-d)-0.0 148)*s,0.0,0.0, 1 .0,0. 0,0.0); 
setClipPlane(-(apw+d)*2*s, -0.0296*5,0.0, 1 .0, 1 .0,0.0); 
setClipPlane(-(apw+d)*2*s,-0.0296*s,0.0,- 1 .0, 1 .0,0.0); 
clipAsAHole(); 

/*Adiskforvisuaiization of blur circle at precision aperture*/ 

makeADisk(0.0,0.0,(2.9404+o+sh)*s,0.0,0.0,-1.0,0.3*s); 

setProperties(1.0,1.0); 

nextOBJ(); 

/*Loop to vary input angles for determination of OPSF*/ 
for (phi=88.6; phi<90.0; phi=phi+0.1) 

{ 

for(theta=0.0; theta < 4.5; theta=theta+0.1) 

{ 

setSourceDir(sin(theta* PI/ 1 80.0)*sin(phi* PI/ 1 80.0),cos(phi*PI/l 80.0), 
cos(theta*PI/180.0)*sin(phi* PI/ 180.0)); 
startRay T race( 1 00000) ; 

} 

} 

} 
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The following code was written in FORTRAN. It is used to post process the output files generated from 
running the provided C code (the first eleven (11) pages of Appendix C). This particular code is used for 
the case in which five detectors, thus five precision apertures are placed at the focal plane of the telescope. 
The output files generated from this code are in a form that can be directly imported into excel, and plotted 
to obtain the Optical Point Spread Function at each detector. Note that all dimensions are in inches and 
are scaled by ten (10). 


PROGRAM OPSF 
IMPLICIT NONE 

REAL*8 APW,D,XP0S,YP0S,ZP0S,AP1,AP2,AP3,MAX1,MAX2,MAX3,MAXNUM1, 
$MAXNUM2,MAXNUM3 

INTEGER I,NUMl,NUM2,NUM3,LINES,J,NUMTHETA,NUMPHI,N,NUMTOTAL,COUNT 
$,NUM,L,P,R,0,C 
CHARACTER* 100 FICHIER 
CHARACTER* 100 SC_LET 

DIMENSION APl(15,45), AP2(15,45), AP3(15,45),MAX1(15), 

$MAX2(15),MAX3(15) 


* Open files for output (apl, ap2, and ap3). Open file numberb.dat which contains the ends 

* of the file names produced by Felix’s code. Numberb.dat contains lines ‘001. dat’, ‘002.dat’, 

* etc. up to the number of output files produced. 

OPEN(15,FILE- numberb.dat', STATUS- OLD') 

OPEN(16,FILE='ap 1 .dat\STATUS='OLD') 

OPEN( 1 7, FILE- ap2.dat', STATUS-OLD') 

OPEN( 1 8,FILE=’ap3 .dat\STATUS=’OLD') 

* APW is the width of the precision apertures, D is the distance between them. 

APW=0.296D0 

D=0. 1 DO 

NUM1=0 

NUM2=0 

NUM3=0 

N=1 

COUNT= 1 

NUMTHETA=45 

NUMPHI=15 

NUMTOTAL=NUMTHETA*NUMPHI 

* Initialize the value in the arrays containing maximum response for each valueof phi (max 1 ...), 

* initialize the maximum response of each aperture (maxnum 1,...) to zero. 

DO 1 C=1,NUMPHI 
MAX1(C)=0 
MAX2(C)=0 
MAX3(C)=0 
1 CONTINUE 
MAXNUM 1=0 
MAXNUM2=0 
MAXNUM3=0 

DO WHILE (N.LE.NUMPHI) 


DO 10 J=COUNT,NUMTHETA*N 
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NUM1=0 

NUM2=0 

NUM3=0 

READ(15,*)SC_LET 

FICHIER = 'OUTPUT 
$ //SC LET 

OPEN(24,FILE=FICHIER,STATUS='OLD') 

* Start a loop in which each output file is opened and read. Each output file corresponds to a given 

* theta, phi combination (see previous C code). When reading a given output file, each 

* line contains a set of coordinates for a ray arriving at the focal plane. If the arriving ray falls within 

* the central aperture, the counter NUM1 is incremented. Similarly, if it falls in the second 

* aperture, NUM2 is incremented and if it falls within the third aperture, NUM3 is incremented. 


DO 235 1=1,50000 

READ(24,*, END=30)XPOS,YPOS,ZPOS 

IF(XPOS.GT.(-APW/2.0DO).AND.XPOS.LE.O)THEN 

!F(ABS(YPOS).LT.(XPOS+0.296dO))THEN 

NUM1=NUM1 + 1 

END IF 

END IF 

IF(XPOS.GT.O.AND.XPOS.LT.(APW/2.dO))THEN 

!F(ABS(YPOS).LT.(0.296dO-XPOS))THEN 

NUM1=NUM1+1 

END IF 

END IF 

IF(ABS(XPOS).GT.(APW/2.dO+D).AND.ABS(XPOS) 

&.LE.(APW+D))THEN 

IF(ABS(YPOS).LT.((abs(XPOS)-(APW+D))+0.296dO))THEN 
NUM2=NUM2+ 1 
END IF 
END IF 

IF(ABS(XPOS).GT.(APW+D).AND.ABS(XPOS).LT. 
&(1.5D0*APW+D))THEN 
IF(ABS(YPOS).LT.(0.296dO-(ABS(XPOS)- 
&(APW+D))))THEN 
NUM2=NUM2+ 1 
END IF 
END IF 

IF(ABS(XPOS).GT.(1.5dO*APW+2.dO*D).AND.ABS(XPOS) 

&.LE.(2.0dO*(APW+D)))THEN 

IF(ABS(YPOS).LT.(ABS(XPOS)-(2.0d0*(APW+D))+0.296d0)) 

&THEN 

NUM3=NUM3+1 
END IF 
END IF 

IF(ABS(XPOS).GT.(2.0dO*(APW+D)).AND.ABS(XPOS).LT. 
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&(2.5dO*APW+2.dO*D))THEN 
IF(ABS(YPOS).LT.(0.296dO-(ABS(XPOS)-(2,OdO* 
&(APW+D)))))THEN 
NUM3=NUM3+ 1 
END IF 
END IF 

NUM= J-NUMTHET A *(N- 1 ) 


235 CONTINUE 

30 AP 1 (N,NUM)=NUM 1 
AP2(N,NUM)=NUM2 
AP3(N,NUM)=NUM3 

10 CONTINUE 


COUNT=J 
N=N+1 
END DO 


* Determine the maximum response within each of the detectors 
DO 60 R= l ,NUMPHI 


70 

60 


DO 70 P= 1 .NUMTHET A 

IF (AP 1 (R,P).GT.MAX 1 (R)) THEN 

MAX 1 (R)=AP 1 (R,P) 

END IF 

IF (AP2(R,P).GT.MAX2(R))THEN 
MAX2(R)=AP2(R,P) 

END IF 

IF (AP3(R,P).GT.MAX3(R))THEN 
MAX3(r)=AP3(R,P) 

END IF 
CONTINUE 
CONTINUE 


DO 80 0=l,NUMPHI 

IF (MAX 1 (O).GT.MAXNUM 1 ) THEN 

MAXNUM 1 =MAX 1(0) 

END IF 

IF (MAX2(0).GT.MAXNUM2) THEN 
MAXNUM2=MAX2(0) 

END IF 

IF (MAX3(0).GT.MAXNUM3) THEN 
MAXNUM3=MAX3(0) 

END IF 

80 CONTINUE 

* Write the output for the Optical Point Spread Function for each detector, normalized 

* by the maximum response within all detectors (maximum response will be at central 

* detector. 

DO 50 L=l, NUMTHET A 

WRITE( 1 6,200) AP 1 ( 1 5,L)/maxnum 1 , AP 1 ( 1 4,L)/maxnum 1 , 



